Significance: Headwater streams account for a majority of river networks worldwide and have a disproportionately large influence on the functioning of aquatic ecosystems. Headwater streams also support critical habitat for many species, including cold-water fishes, many of which are declining or at risk of extinction.
Problem: Headwater streams are largely underrepresented in streamflow monitoring networks, which place greater emphasis on mainstem rivers. As a result, less is known about how headwaters respond to changing water availability. Headwater streams therefore represent a blind spot in understanding flow regime variability and assessing the vulnerability of cold-water fish to changing climatic conditions, including the increasing frequency and severity of drought.
Existential question: How do streamflow regimes vary spatially in headwater stream networks and what does this imply about the sensitivity of these systems (and the species they support) to changing climatic conditions (i.e., drought)?
Objectives:
Conceptual Model: develop a conceptual framework and set of hypotheses based on the objectives above How do we link/combine what we have learned so far to enhance our understanding of spatial diversity in headwater streamflow regimes?
Boxes: Use boxes to highlight additional details of the data, vignettes, and case studies that demonstrate spatiotemporal streamflow heterogeneity
Site information and map
# site information
siteinfo <- read_csv("C:/Users/jbaldock/OneDrive - DOI/Documents/USGS/EcoDrought/EcoDrought Working/Data/EcoDrought_SiteInformation.csv")
siteinfo_sp <- st_as_sf(siteinfo, coords = c("long", "lat"), crs = 4326)
mapview(siteinfo_sp, zcol = "designation")Little g daily data
# flow/yield (and temp) data
dat <- read_csv("C:/Users/jbaldock/OneDrive - DOI/Documents/USGS/EcoDrought/EcoDrought Working/Data/EcoDrought_FlowTempData_DailyWeekly.csv") %>%
mutate(site_name = dplyr::recode(site_name, "Leidy Creek Mouth NWIS" = "Leidy Creek Mouth", "SF Spread Creek Lower NWIS" = "SF Spread Creek Lower", "Dugout Creek NWIS" = "Dugout Creek", "Shields River ab Smith NWIS" = "Shields River Valley Ranch")) %>%
filter(!site_name %in% c("Avery Brook NWIS", "West Brook 0", "BigCreekMiddle", # drop co-located sites
"South River Conway NWIS", "North Fork Flathead River NWIS", # drop big Gs
"Pacific Creek at Moran NWIS", "Shields River nr Livingston NWIS", # drop big Gs
"Donner Blitzen River nr Frenchglen NWIS", # drop big Gs
"WoundedBuckCreek")) %>% # drop little g outside of focal basin
group_by(site_name, basin, subbasin, region, date) %>%
summarize(flow_mean = mean(flow_mean),
tempc_mean = mean(tempc_mean),
Yield_mm = mean(Yield_mm),
Yield_filled_mm = mean(Yield_filled_mm)) %>%
ungroup()
# add water/climate year variables and fill missing dates
dat <- fill_missing_dates(dat, dates = date, groups = site_name)
dat <- add_date_variables(dat, dates = date, water_year_start = 10)Clean and bind little g data (for each basin, restrict to time period for which data quality/availability is ~consistent)
dat_clean <- bind_rows(
dat %>% filter(site_name %in% unlist(siteinfo %>% filter(subbasin == "West Brook") %>% select(site_name)), year(date) >= 2020, date <= date("2025-01-03")) %>%
mutate(Yield_filled_mm = ifelse(site_name == "West Brook Upper" & date > date("2024-10-06"), NA, Yield_filled_mm)) %>%
mutate(Yield_filled_mm = ifelse(site_name == "Mitchell Brook" & date > date("2021-02-28") & date < date("2021-03-26"), NA, Yield_filled_mm)) %>%
mutate(Yield_filled_mm = ifelse(site_name == "Mitchell Brook" & date > date("2021-11-01") & date < date("2022-05-01"), NA, Yield_filled_mm)) %>%
mutate(Yield_filled_mm = ifelse(site_name == "Jimmy Brook" & date > date("2024-12-10"), NA, Yield_filled_mm)) %>%
mutate(subbasin = "West Brook"),
dat %>% filter(site_name %in% unlist(siteinfo %>% filter(subbasin == "Paine Run") %>% select(site_name)), date >= as_date("2018-11-07"), date <= as_date("2023-05-15")) %>% mutate(subbasin = "Paine Run"),
dat %>% filter(site_name %in% unlist(siteinfo %>% filter(subbasin == "Staunton River") %>% select(site_name)), date >= as_date("2018-11-07"), date <= as_date("2022-10-19")) %>% mutate(subbasin = "Staunton River"),
dat %>% filter(site_name %in% c(unlist(siteinfo %>% filter(subbasin == "Big Creek") %>% select(site_name)), "North Fork Flathead River NWIS"), date >= date("2018-08-08"), date <= date("2023-08-03"), site_name != "SkookoleelCreek", Yield_filled_mm > 0) %>% mutate(subbasin = "Big Creek"),
dat %>% filter(site_name %in% c(unlist(siteinfo %>% filter(subbasin == "Coal Creek") %>% select(site_name)), "North Fork Flathead River NWIS"), date >= date("2018-07-29"), date <= date("2023-08-03")) %>% mutate(subbasin = "Coal Creek"),
dat %>% filter(site_name %in% c(unlist(siteinfo %>% filter(subbasin == "McGee Creek") %>% select(site_name)), "North Fork Flathead River NWIS"), date >= date("2017-07-30"), date <= date("2023-12-11")) %>% mutate(subbasin = "McGee Creek"),
dat %>% filter(subbasin == "Snake River", date >= date("2016-04-01"), date <= date("2023-10-03"), site_name != "Leidy Creek Upper") %>% mutate(subbasin = "Snake River"),
dat %>% filter(subbasin == "Shields River", date >= date("2016-04-01"), date <= date("2023-12-31"), site_name != "Brackett Creek") %>%
mutate(logYield = log10(Yield_filled_mm)) %>% mutate(subbasin = "Shields River"),
dat %>% filter(subbasin == "Duck Creek", date >= date("2015-04-01"), date <= date("2023-12-31")) %>% mutate(subbasin = "Duck Creek"),
dat %>% filter(subbasin == "Donner Blitzen", date >= as_date("2019-04-23"), date <= as_date("2022-12-31"), !site_name %in% c("Indian Creek NWIS", "Little Blizten River NWIS")) %>% mutate(subbasin = "Donner Blitzen")
) %>%
filter(Yield_filled_mm > 0) %>%
mutate(logYield = log10(Yield_filled_mm),
designation = "little",
doy_calendar = yday(date)) %>%
select(-Yield_mm) %>%
rename(Yield_mm = Yield_filled_mm)
head(dat_clean)# A tibble: 6 × 16
site_name basin subbasin region date flow_mean tempc_mean Yield_mm
<chr> <chr> <chr> <chr> <date> <dbl> <dbl> <dbl>
1 Avery Brook West Bro… West Br… Mass 2020-01-08 5.96 0.594 1.99
2 Avery Brook West Bro… West Br… Mass 2020-01-09 4.81 0.0336 1.61
3 Avery Brook West Bro… West Br… Mass 2020-01-10 4.88 0.363 1.63
4 Avery Brook West Bro… West Br… Mass 2020-01-11 6.43 1.77 2.15
5 Avery Brook West Bro… West Br… Mass 2020-01-12 21.2 2.81 7.08
6 Avery Brook West Bro… West Br… Mass 2020-01-13 14.3 1.92 4.78
# ℹ 8 more variables: CalendarYear <dbl>, Month <dbl>, MonthName <fct>,
# WaterYear <dbl>, DayofYear <dbl>, logYield <dbl>, designation <chr>,
# doy_calendar <dbl>Load big/super G data
nwis_daily <- read_csv("C:/Users/jbaldock/OneDrive - DOI/Documents/USGS/EcoDrought/EcoDrought Working/Data/EcoDrought_NWIS_FlowTempData_Raw_Daily.csv") %>%
filter(designation == "big",
year(date) >= 1970,
site_name != "Shields River nr Livingston NWIS") %>%
mutate(flowcfs = ifelse(site_name == "Rapidan River NWIS" & date > date("1995-06-26") & date < date("1995-07-01"), NA, flowcfs),
flow_mean_cms = flowcfs*0.02831683199881,
area_sqkm = area_sqmi*2.58999)
# sites
sites <- unique(nwis_daily$site_name)
# site-specific basin area in square km
basinarea <- nwis_daily %>% filter(!is.na(site_id)) %>% group_by(site_name) %>% summarize(area_sqkm = unique(area_sqkm))
# calculate yield
yield_list <- list()
for (i in 1:length(sites)) {
d <- nwis_daily %>% filter(site_name == sites[i])
ba <- unlist(basinarea %>% filter(site_name == sites[i]) %>% select(area_sqkm))
yield_list[[i]] <- add_daily_yield(data = d, values = flow_mean_cms, basin_area = as.numeric(ba))
}
nwis_daily_wyield <- do.call(rbind, yield_list)
dat_clean_big <- nwis_daily_wyield %>%
select(site_name, basin, subbasin, region, date, Yield_mm, tempc, flowcfs) %>%
mutate(logYield = log10(Yield_mm), doy_calendar = yday(date)) %>%
rename(tempc_mean = tempc, flow_mean = flowcfs)
# add water/climate year variables and fill missing dates
dat_clean_big <- fill_missing_dates(dat_clean_big, dates = date, groups = site_name)
dat_clean_big <- add_date_variables(dat_clean_big, dates = date, water_year_start = 10)
head(dat_clean_big)# A tibble: 6 × 15
site_name basin subbasin region date Yield_mm tempc_mean flow_mean
<chr> <chr> <chr> <chr> <date> <dbl> <dbl> <dbl>
1 South River Co… West… West Br… Mass 1970-01-01 1.81 NA 46
2 South River Co… West… West Br… Mass 1970-01-02 1.69 NA 43
3 South River Co… West… West Br… Mass 1970-01-03 1.61 NA 41
4 South River Co… West… West Br… Mass 1970-01-04 1.53 NA 39
5 South River Co… West… West Br… Mass 1970-01-05 1.49 NA 38
6 South River Co… West… West Br… Mass 1970-01-06 1.49 NA 38
# ℹ 7 more variables: logYield <dbl>, doy_calendar <dbl>, CalendarYear <dbl>,
# Month <dbl>, MonthName <fct>, WaterYear <dbl>, DayofYear <dbl>View big G streamflow time series data
dat_clean_big %>% ggplot() + geom_line(aes(x = date, y = logYield)) + facet_wrap(~site_name, nrow = 8) + theme_bw()
Download Daymet precip data and summarize by water year
# big G site lat/long
mysites <- nwis_daily %>% group_by(site_name, basin, subbasin, region) %>% summarize(lat = unique(lat), long = unique(long)) %>% ungroup()
# download point location Daymet data
climlist <- vector("list", length = dim(mysites)[1])
for (i in 1:dim(mysites)[1]) {
clim <- download_daymet(site = mysites$site_name[i], lat = mysites$lat[i], lon = mysites$long[i], start = 1980, end = 2024, internal = T)
climlist[[i]] <- tibble(clim$data) %>%
mutate(air_temp_mean = (tmax..deg.c. + tmin..deg.c.)/2,
date = as.Date(paste(year, yday, sep = "-"), "%Y-%j"),
site_name = mysites$site_name[i]) %>%
select(12,2,11,10,4,6) %>% rename(precip_mmday = 5, swe_kgm2 = 6)
print(i)
}[1] 1[1] 2[1] 3[1] 4[1] 5[1] 6[1] 7[1] 8# combine and add water years
climdf <- do.call(rbind, climlist) %>% left_join(mysites) %>% mutate(year = year(date))
climdf <- add_date_variables(climdf, dates = date, water_year_start = 10)
# calculate total annual precipitation in mm, by site and water year
climdf_summ <- climdf %>%
group_by(site_name, basin, subbasin, region, WaterYear) %>%
summarize(precip_total = sum(precip_mmday), sampsize = n()) %>%
mutate(precip_total_z = scale(precip_total)) %>%
ungroup() %>%
filter(sampsize >= 350)View time series and scatter plots of big G water availability (total annual yield, sum of daily values, black lines) and precipitation (total annual precip, blues lines). Note that the panel labels indicate basins, not NWIS gage names. In the manuscript, pair these plots with CONUS map and detailed inset maps of focal basins (Figure 1)
# calculate annual water availability at big g as sum of daily yield values...retain only years with >95% data coverage (at least 350 days)
wateravail <- dat_clean_big %>%
filter(!is.na(Yield_mm)) %>%
group_by(basin, site_name, WaterYear) %>%
summarize(sampsize = n(), totalyield = sum(Yield_mm, na.rm = TRUE)) %>%
mutate(totalyield_z = scale(totalyield)[,1]) %>%
ungroup() %>%
filter(sampsize >= 350) %>%
complete(basin, WaterYear = 1971:2024, fill = list(sampsize = NA, totalyield = NA))
# get range of years for little g data
daterange <- dat_clean %>% group_by(basin) %>% summarize(minyear = year(min(date)), maxyear = year(max(date)))
# spread ecod years
mylist <- vector("list", length = dim(daterange)[1])
for (i in 1:dim(daterange)[1]) {
mylist[[i]] <- tibble(basin = daterange$basin[i], WaterYear = seq(from = daterange$minyear[i], to = daterange$maxyear[i], by = 1))
}
yrdf <- do.call(rbind, mylist) %>% mutate(ecodyr = "yes")ggplot() +
geom_rect(data = daterange, aes(xmin = minyear-0.5, xmax = maxyear+0.5, ymin = -Inf, ymax = +Inf), fill = "grey") +
geom_line(data = wateravail, aes(x = WaterYear, y = totalyield), linewidth = 1) +
geom_line(data = climdf_summ, aes(x = WaterYear, y = precip_total), linewidth = 0.5, col = "blue") +
facet_wrap(~basin) +
xlab("Water year") + ylab("Total annual yield (mm) / Total annual precipitation (mm)") +
theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggplot() +
geom_rect(data = daterange, aes(xmin = minyear-0.5, xmax = maxyear+0.5, ymin = -Inf, ymax = +Inf), fill = "grey") +
geom_line(data = wateravail, aes(x = WaterYear, y = totalyield_z), linewidth = 1) +
geom_line(data = climdf_summ, aes(x = WaterYear, y = precip_total_z), linewidth = 0.5, col = "blue") +
facet_wrap(~basin) +
xlab("Water year") + ylab("Total annual yield (scaled) / Total annual precipitation (scaled)") +
theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
wateravail %>% select(-sampsize) %>% left_join(climdf_summ %>% select(-sampsize)) %>% left_join(yrdf) %>%
ggplot(aes(x = precip_total_z, y = totalyield_z)) +
geom_point(aes(color = ecodyr)) +
facet_wrap(~basin) +
xlab("Total annual precipitation (scaled)") + ylab("Total annual yield (scaled)") +
theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") +
stat_cor(method = "pearson", aes(label = ..r.label..))
Objective 1: Evaluate the extent and magnitude of spatial and temporal variation in headwater streamflow.
Approach:
Notes:
View daily time series data by sub-basin. Note that we are using the “Super G” NWIS data for the reference gage (black line). Per Robert comment, entirely nested design is cute, but doesn’t reflect how the data is actually used.
Big G NWIS sites/reference gages for each basin/subbasin:
dat_clean_big %>% group_by(region, basin, subbasin) %>% summarize(site_name = unique(site_name)) %>% ungroup() %>% filter(!is.na(region)) %>% kable()| region | basin | subbasin | site_name |
|---|---|---|---|
| Flat | Flathead | Flathead | North Fork Flathead River NWIS |
| Mass | West Brook | West Brook | South River Conway NWIS |
| Oreg | Donner Blitzen | Donner Blitzen | Donner Blitzen River nr Frenchglen NWIS |
| Shen | Paine Run | Paine Run | South River Harriston NWIS |
| Shen | Piney River | Piney River | Battle Run NWIS |
| Shen | Staunton River | Staunton River | Rapidan River NWIS |
| Shields | Shields River | Shields River | Yellowstone River Livingston NWIS |
| Snake | Snake River | Snake River | Pacific Creek at Moran NWIS |
Plotting function
# set up color palette
mycols <- c(brewer.pal(8, "Dark2"), "dodgerblue", "darkorchid")
# create plotting function
myplotfun <- function(subbas, lab, bigG) {
data_sub <- dat_clean %>% filter(subbasin == subbas)
par(mar = c(1.7,1.7,0.1,0.1), mgp = c(2.5,0.6,0))
plot(logYield ~ date, data_sub, type = "n", xlab = "", ylab = "")
#tempdat_little <- data_sub %>% filter(designation == "little")
tempdat_little <- fill_missing_dates(data_sub, dates = date, groups = site_name, pad_ends = FALSE)
tempdat_big <- dat_clean_big %>% filter(site_name == bigG, date >= min(tempdat_little$date), date <= max(tempdat_little$date))
#tempdat_big <- data_sub %>% filter(designation == "big")
tempdat_big <- fill_missing_dates(tempdat_big, dates = date, groups = site_name, pad_ends = FALSE)
mysites <- sort(unique(tempdat_little$site_name))
for (i in 1:length(mysites)) {
lines(logYield ~ date, tempdat_little %>% filter(site_name == mysites[i]), col = mycols[i], lwd = 0.7)
}
lines(logYield ~ date, tempdat_big, col = "white", lwd = 1.7)
lines(logYield ~ date, tempdat_big, col = "black", lwd = 1)
usr <- par("usr")
par(usr = c(0,1,0,1))
text(0.02, 0.9, labels = lab, cex = 1.2)
par(usr = usr)
}Generate plot
jpeg("C:/Users/jbaldock/OneDrive - DOI/Documents/USGS/EcoDrought/EcoDrought Working/EcoDrought-Analysis/Qualitative/EcoD_timeseries.jpg", width = 8, height = 12, units = "in", res = 1000)
par(mfrow = c(10,1), oma = c(0.5,2,0.5,0.5))
# West Brook
myplotfun(subbas = "West Brook", lab = "(a)", bigG = "South River Conway NWIS")
# Paine
myplotfun(subbas = "Paine Run", lab = "(b)", bigG = "South River Harriston NWIS")
# Staunton River
myplotfun(subbas = "Staunton River", lab = "(c)", bigG = "Rapidan River NWIS")
# Big Creek
myplotfun(subbas = "Big Creek", lab = "(d)", bigG = "North Fork Flathead River NWIS")
# Coal Creek
myplotfun(subbas = "Coal Creek", lab = "(e)", bigG = "North Fork Flathead River NWIS")
# McGee Creek
myplotfun(subbas = "McGee Creek", lab = "(f)", bigG = "North Fork Flathead River NWIS")
# Snake River
myplotfun(subbas = "Snake River", lab = "(g)", bigG = "Pacific Creek at Moran NWIS")
# Shields River
myplotfun(subbas = "Shields River", lab = "(h)", bigG = "Yellowstone River Livingston NWIS")
# Duck Creek
myplotfun(subbas = "Duck Creek", lab = "(i)", bigG = "Yellowstone River Livingston NWIS")
# Donner Blitzen
myplotfun(subbas = "Donner Blitzen", lab = "(j)", bigG = "Donner Blitzen River nr Frenchglen NWIS")
# common axis label
mtext("log(Daily yield, mm)", outer = TRUE, side = 2, line = 0.5)
dev.off()png
2 
these all should probably be trimmed to just the summer data…or comparable period over which data availability is maximized across sites and subbasins
Need to do
Plotting functions
# create plotting function
myplotfun_ex <- function(subbas, lab) {
# filter to subbasin and get site designation
data_sub <- dat_clean %>% filter(subbasin == subbas, Month %in% c(7,8,9))
sitesdesig <- data_sub %>% group_by(site_name) %>% summarize(designation = unique(designation)) %>% ungroup()
# calculate exceedance probability by site
exeeddat <- data_sub %>%
select(date, site_name, logYield) %>%
spread(key = site_name, value = logYield) %>%
drop_na() %>%
gather(key = site_name, value = logYield, 2:ncol(.)) %>%
group_by(site_name) %>%
arrange(desc(logYield), .by_group = TRUE) %>%
mutate(exceedance = 100/length(logYield)*1:length(logYield)) %>%
ungroup() %>%
left_join(sitesdesig)
exeeddat_little <- exeeddat %>% filter(designation == "little")
exeeddat_big <- exeeddat %>% filter(designation == "big")
# set up plot
par(mar = c(1.7,1.7,0.1,0.1), mgp = c(2.5,0.6,0))
mysites <- sort(unique(exeeddat_little$site_name))
plot(logYield ~ exceedance, exeeddat, type = "n", xlab = "", ylab = "")
# add little g's
for (i in 1:length(mysites)) {
lines(logYield ~ exceedance, exeeddat_little %>% filter(site_name == mysites[i]), col = mycols[i], lwd = 1)
}
# add big G
lines(logYield ~ exceedance, exeeddat_big, col = "white", lwd = 2.7)
lines(logYield ~ exceedance, exeeddat_big, col = "black", lwd = 2)
# panel label
usr <- par("usr")
par(usr = c(0,1,0,1))
text(0.92, 0.9, labels = lab, cex = 1.2)
par(usr = usr)
}
# CV of little g Yield by exceedance probability
myplotfun_cv <- function(subbas, lab) {
# filter to subbasin and get site designation
data_sub <- dat_clean %>% filter(subbasin == subbas)
sitesdesig <- data_sub %>% group_by(site_name) %>% summarize(designation = unique(designation)) %>% ungroup()
# calculate exceedance probability by site
exeeddat <- data_sub %>%
select(date, site_name, logYield) %>%
spread(key = site_name, value = logYield) %>%
drop_na() %>%
gather(key = site_name, value = logYield, 2:ncol(.)) %>%
group_by(site_name) %>%
arrange(desc(logYield), .by_group = TRUE) %>%
mutate(exceedance = 100/length(logYield)*1:length(logYield)) %>%
ungroup() %>%
left_join(sitesdesig)
exeeddat_little <- exeeddat %>% filter(designation == "little")
exeeddat_big <- exeeddat %>% filter(designation == "big")
# set up plot
par(mar = c(1.7,1.7,0.1,0.1), mgp = c(2.5,0.6,0))
mysites <- sort(unique(exeeddat_little$site_name))
tiblist <- list()
for (i in 1:length(mysites)) {
tt <- exeeddat %>% filter(site_name == mysites[i])
mylinint <- approx(x = tt$exceedance, y = tt$logYield, xout = seq(from = 0, to = 100, by = 1))
tiblist[[i]] <- tibble(site_name = mysites[i], exceedance = mylinint$x, logYield = mylinint$y)
}
cvtib <- do.call(bind_rows, tiblist) %>%
group_by(exceedance) %>%
summarize(sdf = sd(logYield)) %>%
ungroup()
plot(sdf ~ exceedance, cvtib, type = "l", col = "grey40", lwd = 2, ylim = c(0,0.8))
# panel label
usr <- par("usr")
par(usr = c(0,1,0,1))
text(0.92, 0.9, labels = lab, cex = 1.2)
par(usr = usr)
} Plot exceedance
jpeg("C:/Users/jbaldock/OneDrive - DOI/Documents/USGS/EcoDrought/EcoDrought Working/EcoDrought-Analysis/Qualitative/EcoD_exceedance_JAS.jpg", width = 4, height = 8, units = "in", res = 1000)
par(mfrow = c(5,2), oma = c(2,2,0.5,0.5))
myplotfun_ex(subbas = "West Brook", lab = "(a)")
myplotfun_ex(subbas = "Paine Run", lab = "(b)")
myplotfun_ex(subbas = "Staunton River", lab = "(c)")
myplotfun_ex(subbas = "Big Creek", lab = "(d)")
myplotfun_ex(subbas = "Coal Creek", lab = "(e)")
myplotfun_ex(subbas = "McGee Creek", lab = "(f)")
myplotfun_ex(subbas = "Snake River", lab = "(g)")
myplotfun_ex(subbas = "Shields River", lab = "(h)")
myplotfun_ex(subbas = "Duck Creek", lab = "(i)")
myplotfun_ex(subbas = "Donner Blitzen", lab = "(j)")
mtext("log(Daily yield, mm)", outer = TRUE, side = 2, line = 0.5)
mtext("Exceedance probability", outer = TRUE, side = 1, line = 0.5)
dev.off()
#Plot CV by exceedance
jpeg("C:/Users/jbaldock/OneDrive - DOI/Documents/USGS/EcoDrought/EcoDrought Working/EcoDrought-Analysis/Qualitative/EcoD_exceedance_cv.jpg", width = 4, height = 8, units = "in", res = 1000)
par(mfrow = c(5,2), oma = c(2,2,0.5,0.5))
myplotfun_cv(subbas = "West Brook", lab = "(a)")
myplotfun_cv(subbas = "Paine Run", lab = "(b)")
myplotfun_cv(subbas = "Staunton River", lab = "(c)")
myplotfun_cv(subbas = "Big Creek", lab = "(d)")
myplotfun_cv(subbas = "Coal Creek", lab = "(e)")
myplotfun_cv(subbas = "McGee Creek", lab = "(f)")
myplotfun_cv(subbas = "Snake River", lab = "(g)")
myplotfun_cv(subbas = "Shields River", lab = "(h)")
myplotfun_cv(subbas = "Duck Creek", lab = "(i)")
myplotfun_cv(subbas = "Donner Blitzen", lab = "(j)")
mtext("SD of log(daily yield, mm) at little g's", outer = TRUE, side = 2, line = 0.5)
mtext("Exceedance probability", outer = TRUE, side = 1, line = 0.5)
dev.off()
#Under the current objectives/framework, I’m not sure there is much of a role for these plots. Leaving in for reference
myridgesfun <- function(subbas, bigG) {
td <- dat_clean %>% filter(subbasin == subbas)
td2 <- td %>%
group_by(subbasin, site_name, designation, CalendarYear, Month, MonthName) %>%
summarize(logYield = mean(logYield)) %>%
ungroup() %>%
mutate(MonthName = factor(MonthName, levels = rev(c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))),
CalendarYear = factor(CalendarYear))
tempdat_big <- dat_clean_big %>%
filter(site_name == bigG, date >= min(td$date), date <= max(td$date)) %>%
group_by(site_name, CalendarYear, Month, MonthName) %>%
summarize(logYield = mean(logYield)) %>%
ungroup() %>%
mutate(MonthName = factor(MonthName, levels = rev(c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))), CalendarYear = factor(CalendarYear))
return(ggplot(data = td2) +
geom_density_ridges(data = td2, aes(x = logYield, y = MonthName), alpha = 0.5,point_alpha = 0.2) +
geom_line(data = td2, aes(x = logYield, y = MonthName, group = site_name, color = site_name), orientation = "y", alpha = 0.3) +
geom_point(data = td2, aes(x = logYield, y = MonthName, group = site_name, color = site_name), alpha = 0.3) +
geom_line(data = tempdat_big, aes(x = logYield, y = MonthName, group = site_name), orientation = "y", alpha = 0.5) +
geom_point(data =tempdat_big, aes(x = logYield, y = MonthName, group = site_name), alpha = 0.5) +
theme_bw() + theme(legend.position = "none") +
facet_wrap2(~CalendarYear, nrow = 3, ncol = 3, trim_blank = FALSE) +
xlab("Monthly mean log(Yield, mm/day)") + ylab(""))
}myridgesfun(subbas = "West Brook", bigG = "South River Conway NWIS")
myridgesfun(subbas = "Paine Run", bigG = "South River Harriston NWIS")
myridgesfun(subbas = "Staunton River", bigG = "Rapidan River NWIS")
myridgesfun(subbas = "Big Creek", bigG = "North Fork Flathead River NWIS")
myridgesfun(subbas = "Coal Creek", bigG = "North Fork Flathead River NWIS")
myridgesfun(subbas = "McGee Creek", bigG = "North Fork Flathead River NWIS")
myridgesfun(subbas = "Snake River", bigG = "Pacific Creek at Moran NWIS")
myridgesfun(subbas = "Shields River", bigG = "Yellowstone River Livingston NWIS")
myridgesfun(subbas = "Duck Creek", bigG = "Yellowstone River Livingston NWIS")
myridgesfun(subbas = "Donner Blitzen", bigG = "Donner Blitzen River nr Frenchglen NWIS")
Purpose: Explore non-linear and hysteretic relationship between flow at reference vs. at headwater gages
Approach:
Find little g site/water years with at least 90% data completeness, join big and little g daily yield data
# find little g site/water yeras with at least 90% data availability
mysiteyrs <- dat_clean %>%
group_by(site_name, basin, subbasin, region, WaterYear) %>%
summarize(numdays = n(),
totalyield_site = sum(Yield_mm)) %>%
ungroup() %>%
filter(numdays >= 0.9*365) %>%
mutate(forhyst = 1) %>%
left_join(wateravail %>% select(basin, WaterYear, totalyield_z))
# join big G data to little g
dat_clean_hyst <- dat_clean %>%
left_join(mysiteyrs) %>%
filter(forhyst == 1) %>%
left_join(dat_clean_big %>% select(basin, date, Yield_mm, logYield) %>% rename(Yield_mm_big = Yield_mm, logYield_big = logYield)) %>%
left_join(wateravail %>% select(basin, WaterYear, totalyield_z))
# view data
mysiteyrs# A tibble: 87 × 9
site_name basin subbasin region WaterYear numdays totalyield_site forhyst
<chr> <chr> <chr> <chr> <dbl> <int> <dbl> <dbl>
1 Avery Brook West… West Br… Mass 2021 365 872. 1
2 Avery Brook West… West Br… Mass 2022 365 602. 1
3 Avery Brook West… West Br… Mass 2023 365 1147. 1
4 Avery Brook West… West Br… Mass 2024 345 993. 1
5 Big Creek NW… Flat… Big Cre… Flat 2019 329 566. 1
6 Big Creek NW… Flat… Big Cre… Flat 2020 347 719. 1
7 Big Creek NW… Flat… Big Cre… Flat 2021 352 554. 1
8 Big Creek NW… Flat… Big Cre… Flat 2022 340 842. 1
9 Buck Creek Shie… Shields… Shiel… 2022 329 174. 1
10 CycloneCreek… Flat… Coal Cr… Flat 2019 354 374. 1
# ℹ 77 more rows
# ℹ 1 more variable: totalyield_z <dbl>dat_clean_hyst# A tibble: 30,889 × 22
site_name basin subbasin region date flow_mean tempc_mean Yield_mm
<chr> <chr> <chr> <chr> <date> <dbl> <dbl> <dbl>
1 Avery Brook West Br… West Br… Mass 2020-10-01 1.78 13.5 0.595
2 Avery Brook West Br… West Br… Mass 2020-10-02 1.02 12.5 0.341
3 Avery Brook West Br… West Br… Mass 2020-10-03 0.863 11.0 0.288
4 Avery Brook West Br… West Br… Mass 2020-10-04 0.654 10.0 0.218
5 Avery Brook West Br… West Br… Mass 2020-10-05 0.613 10.2 0.204
6 Avery Brook West Br… West Br… Mass 2020-10-06 0.628 11.1 0.209
7 Avery Brook West Br… West Br… Mass 2020-10-07 0.717 11.4 0.239
8 Avery Brook West Br… West Br… Mass 2020-10-08 0.787 10.7 0.262
9 Avery Brook West Br… West Br… Mass 2020-10-09 0.661 9.04 0.221
10 Avery Brook West Br… West Br… Mass 2020-10-10 0.734 9.95 0.245
# ℹ 30,879 more rows
# ℹ 14 more variables: CalendarYear <dbl>, Month <dbl>, MonthName <fct>,
# WaterYear <dbl>, DayofYear <dbl>, logYield <dbl>, designation <chr>,
# doy_calendar <dbl>, numdays <int>, totalyield_site <dbl>, forhyst <dbl>,
# totalyield_z <dbl>, Yield_mm_big <dbl>, logYield_big <dbl>Create plotting function
hystplotfun <- function(mysite, wy, months = c(1:12)) {
(dat_clean_hyst %>%
filter(site_name == mysite, WaterYear == wy, Month %in% months) %>%
ggplot(aes(x = logYield_big, y = logYield, color = date)) +
geom_segment(aes(xend = c(tail(logYield_big, n = -1), NA), yend = c(tail(logYield, n = -1), NA)), arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_point() +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
scale_color_gradientn(colors = cet_pal(250, name = "c2"), trans = "date") +
xlim(-1,1.5) + ylim(-1.5,1.8) +
#scale_color_scico(palette = "romaO") +
xlab("log(Yield, mm) at reference gage") + ylab("log(Yield, mm) at headwater gage") +
# geom_text(data = annotations, aes(x = xpos, y = ypos, label = annotateText), hjust = "inward", vjust = "inward") +
annotate(geom = "text", x = -Inf, y = Inf, label = paste(mysite, ", WY ", wy, sep = ""), hjust = -0.1, vjust = 1.5) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "grey85"),
axis.title = element_blank()) +
theme(plot.margin = margin(0.1,0.1,0,0, "cm")) )
}
# change to linear color palette
hystplotfunlin <- function(mysite, wy, months = c(1:12)) {
(dat_clean_hyst %>%
filter(site_name == mysite, WaterYear == wy, Month %in% months) %>%
ggplot(aes(x = logYield_big, y = logYield, color = date)) +
geom_segment(aes(xend = c(tail(logYield_big, n = -1), NA), yend = c(tail(logYield, n = -1), NA)), arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_point() +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
scale_color_gradientn(colors = cet_pal(90, name = "r1"), trans = "date") +
xlim(-1,1.5) + ylim(-1.5,1.8) +
#scale_color_scico(palette = "romaO") +
xlab("log(Yield, mm) at reference gage") + ylab("log(Yield, mm) at headwater gage") +
# geom_text(data = annotations, aes(x = xpos, y = ypos, label = annotateText), hjust = "inward", vjust = "inward") +
annotate(geom = "text", x = -Inf, y = Inf, label = paste(mysite, ", WY ", wy, sep = ""), hjust = -0.1, vjust = 1.5) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "grey85"),
axis.title = element_blank()) +
theme(plot.margin = margin(0.1,0.1,0,0, "cm")) )
}Plot example sites and years. Generally, in rain-dominated basins the East (top two rows), we see relatively little hysteresis/non-stationarity in the relationship between streamflow in headwaters and at reference gages. In contrast, in snowmelt-dominated basins of the Rocky Mountains, we see much stronger hysteresis/non-stationarity in the relationship between streamflow in headwaters and at reference gages, but this varies considerably among locations.
ggarrange(hystplotfun(mysite = "West Brook NWIS", wy = 2021) + theme(legend.position = "none"),
hystplotfun(mysite = "Jimmy Brook", wy = 2021) + theme(legend.position = "none"),
hystplotfun(mysite = "Sanderson Brook", wy = 2021) + theme(legend.position = "none"),
hystplotfun(mysite = "Avery Brook", wy = 2021) + theme(legend.position = "none"),
hystplotfun(mysite = "Staunton River 10", wy = 2021) + theme(legend.position = "none"),
hystplotfun(mysite = "Staunton River 06", wy = 2021) + theme(legend.position = "none"),
hystplotfun(mysite = "Staunton River 03", wy = 2021) + theme(legend.position = "none"),
hystplotfun(mysite = "Staunton River 02", wy = 2021) + theme(legend.position = "none"),
hystplotfun(mysite = "Big Creek NWIS", wy = 2020) + theme(legend.position = "none"),
hystplotfun(mysite = "WernerCreek", wy = 2020) + theme(legend.position = "none"),
hystplotfun(mysite = "CycloneCreekUpper", wy = 2020) + theme(legend.position = "none"),
hystplotfun(mysite = "McGeeCreekTrib", wy = 2020) + theme(legend.position = "none"),
hystplotfun(mysite = "Buck Creek", wy = 2022) + theme(legend.position = "none"),
hystplotfun(mysite = "Dugout Creek", wy = 2022) + theme(legend.position = "none"),
hystplotfun(mysite = "EF Duck Creek be HF", wy = 2022) + theme(legend.position = "none"),
get_legend(hystplotfun(mysite = "EF Duck Creek be HF", wy = 2022)),
nrow = 4, ncol = 4) %>%
annotate_figure(left = text_grob("log(Yield, mm) at headwater gage", rot = 90),
bottom = text_grob("log(Yield, mm) at reference gage"))
What if we just consider the summer (low flow) period, July - September? In some ways, we see an opposite pattern relative to that describe above. In that we see greater hysteretic behavior in rain-dominated basins (driven by frequent summer storms and lagged runoff response) and less hysteretic behavior in snow-dominated basins, as all streams are in recession from the spring snowmelt peak. However, in general, we see greater divergence from a 1:1 relationship in snowmelt vs rain-dominated basins.
ggarrange(hystplotfunlin(mysite = "West Brook NWIS", wy = 2021, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "Jimmy Brook", wy = 2021, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "Sanderson Brook", wy = 2021, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "Avery Brook", wy = 2021, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "Staunton River 10", wy = 2021, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "Staunton River 06", wy = 2021, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "Staunton River 03", wy = 2021, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "Staunton River 02", wy = 2021, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "Big Creek NWIS", wy = 2020, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "WernerCreek", wy = 2020, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "CycloneCreekUpper", wy = 2020, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "McGeeCreekTrib", wy = 2020, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "Buck Creek", wy = 2022, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "Dugout Creek", wy = 2022, months = c(7:9)) + theme(legend.position = "none"),
hystplotfunlin(mysite = "EF Duck Creek be HF", wy = 2022, months = c(7:9)) + theme(legend.position = "none"),
get_legend(hystplotfunlin(mysite = "EF Duck Creek be HF", wy = 2022, months = c(7:9))),
nrow = 4, ncol = 4) %>%
annotate_figure(left = text_grob("log(Yield, mm) at headwater gage", rot = 90),
bottom = text_grob("log(Yield, mm) at reference gage"))
ggplotly(hystplotfun(mysite = "West Brook NWIS", wy = 2021))ggplotly(hystplotfun(mysite = "Jimmy Brook", wy = 2021))ggplotly(hystplotfun(mysite = "Sanderson Brook", wy = 2021))ggplotly(hystplotfun(mysite = "Avery Brook", wy = 2021))ggplotly(hystplotfun(mysite = "Staunton River 10", wy = 2021))ggplotly(hystplotfun(mysite = "Staunton River 06", wy = 2021))ggplotly(hystplotfun(mysite = "Staunton River 03", wy = 2021))ggplotly(hystplotfun(mysite = "Staunton River 02", wy = 2021))ggplotly(hystplotfun(mysite = "Big Creek NWIS", wy = 2020))ggplotly(hystplotfun(mysite = "WernerCreek", wy = 2020))ggplotly(hystplotfun(mysite = "CycloneCreekUpper", wy = 2020))ggplotly(hystplotfun(mysite = "McGeeCreekTrib", wy = 2020))ggplotly(hystplotfun(mysite = "Buck Creek", wy = 2022))ggplotly(hystplotfun(mysite = "Dugout Creek", wy = 2022))ggplotly(hystplotfun(mysite = "EF Duck Creek be HF", wy = 2022))Conduct baseflow extraction, event delineation, and calculate event-specific hysteresis index using dynamic time warping (Xue et al. 2024).
# set baseflow extraction and event delineation parameters
alp <- 0.95
numpass <- 3
thresh <- 0.75
# empty list to story site-year data
hysteresis_list <- list()
for (j in 1:nrow(mysiteyrs)) {
# filter data by site name and water year
dd <- dat_clean_hyst %>% filter(site_name == mysiteyrs$site_name[j], WaterYear == mysiteyrs$WaterYear[j])
# baseflow extraction on big G yield
dd <- dd %>%
filter(!is.na(Yield_mm_big)) %>%
mutate(bf = baseflowB(Yield_mm_big, alpha = alp, passes = numpass)$bf,
bfi = baseflowB(Yield_mm_big, alpha = alp, passes = numpass)$bfi) %>%
ungroup()
# delineate events
events <- eventBaseflow(dd$Yield_mm_big, BFI_Th = thresh, bfi = dd$bfi)
events <- events %>% mutate(len = end - srt + 1)
#print("Done - baseflow extraction and event delineation")
## Tidy events
# define positions of non-events
# srt <- c(1)
# end <- c(events$srt[1]-1)
# for (i in 2:(dim(events)[1])) {
# srt[i] <- events$end[i-1]+1
# end[i] <- events$srt[i]-1
# }
# nonevents <- data.frame(tibble(srt, end) %>%
# mutate(len = end - srt) %>%
# filter(len >= 0) %>% select(-len) %>%
# add_row(srt = events$end[dim(events)[1]]+1, end = dim(dd)[1])
# )
# create vectors of binary event/non-event and event IDs
isevent_vec <- rep(2, times = dim(dd)[1])
eventid_vec <- rep(NA, times = dim(dd)[1])
for (i in 1:dim(events)[1]) {
isevent_vec[c(events[i,1]:events[i,2])] <- 1
eventid_vec[c(events[i,1]:events[i,2])] <- i
}
# # create vector of non-event IDs
# noneventid_vec <- rep(NA, times = dim(dd)[1])
# for (i in 1:dim(nonevents)[1]) { noneventid_vec[c(nonevents[i,1]:nonevents[i,2])] <- i }
#
# # create vector of "agnostic events": combined hydro events and non-events
# agnevents <- rbind(events %>% select(srt, end) %>% mutate(event = 1), nonevents %>% mutate(event = 0)) %>% arrange((srt))
# agneventid_vec <- c()
# for (i in 1:dim(agnevents)[1]){ agneventid_vec[c(agnevents[i,1]:agnevents[i,2])] <- i }
# add event/non-event vectors to Big G data
dd <- dd %>%
mutate(isevent = isevent_vec,
eventid = eventid_vec,
#noneventid = noneventid_vec,
#agneventid = agneventid_vec,
big_event_yield = ifelse(isevent_vec == 1, Yield_mm_big, NA),
big_nonevent_yield = ifelse(isevent_vec == 2, Yield_mm_big, NA),
big_event_quick = big_event_yield - bf) %>%
rename(big_bf = bf, big_bfi = bfi)
#print("Done - tidy events")
# calculate hysteresis index for each event using dynamic time warping
event_list <- list()
for (i in 1:max(dd$eventid, na.rm = TRUE)) {
# filter by each event and normalized big and little Yield
ddd <- dd %>% filter(eventid == i)
dddd <- dd %>% filter(date %in% c(min(ddd$date)-1, max(ddd$date)+1))
ddd <- ddd %>%
bind_rows(dddd) %>%
arrange(date) %>%
mutate(weight = (logYield_big - min(logYield_big)) / (max(logYield_big) - min(logYield_big)),
yield_little_norm = (logYield - min(logYield)) / (max(logYield) - min(logYield)),
yield_big_norm = (logYield_big - min(logYield_big)) / (max(logYield_big) - min(logYield_big)))
# ddd %>%
# ggplot(aes(x = yield_big_norm, y = yield_little_norm, color = date)) +
# geom_segment(aes(xend = c(tail(yield_big_norm, n = -1), NA), yend = c(tail(yield_little_norm, n = -1), NA)),
# arrow = arrow(length = unit(0.2, "cm")), color = "black") +
# geom_point() +
# geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
# scale_color_gradientn(colors = cet_pal(250, name = "r1"), trans = "date") +
# theme_bw() +
# theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# align data using dynamic time warping, allow for error handling (e.g., skips if either big or little g is constant during event, b/c normalization returns NaN)
skip_to_next <- FALSE
tryCatch(align <- dtw(x = unlist(ddd %>% select(yield_big_norm)), y = unlist(ddd %>% select(yield_little_norm)), step = asymmetric, keep = TRUE), error = function(e) { skip_to_next <<- TRUE })
if(skip_to_next) { next }
#align <- dtw(x = unlist(ddd %>% select(yield_big_norm)), y = unlist(ddd %>% select(yield_little_norm)), step = asymmetric, keep = TRUE)
# plot(align, type = "threeway")
# plot(align, type = "twoway", offset = -1)
# find time-warped distance, S_i in Xue et al. (2024)
ddd <- ddd %>% mutate(distance = align$index1 - align$index2)
# ddd %>%
# ggplot(aes(x = date, y = distance*weight)) +
# geom_bar(stat = "identity") +
# theme_bw()
# summarize event
event_list[[i]] <- ddd %>%
group_by(site_name, basin, subbasin, region, WaterYear) %>%
summarize(eventdays = n(),
eventid = i,
mindate = min(date),
maxdate = max(date),
totalyieldevent_little = sum(Yield_mm),
totalyieldevent_big = sum(Yield_mm_big),
hysteresis = sum(ddd$weight * ddd$distance) / sum(ddd$weight) # hysteresis index as in Xue et al (2024)
) %>%
ungroup()
}
#print("Done - dynamic time warping")
hysteresis_list[[j]] <- do.call(bind_rows, event_list)
print(j)
}[1] 1[1] 2[1] 3[1] 4[1] 5[1] 6[1] 7[1] 8[1] 9[1] 10[1] 11[1] 12[1] 13[1] 14[1] 15[1] 16[1] 17[1] 18[1] 19[1] 20[1] 21[1] 22[1] 23[1] 24[1] 25[1] 26[1] 27[1] 28[1] 29[1] 30[1] 31[1] 32[1] 33[1] 34[1] 35[1] 36[1] 37[1] 38[1] 39[1] 40[1] 41[1] 42[1] 43[1] 44[1] 45[1] 46[1] 47[1] 48[1] 49[1] 50[1] 51[1] 52[1] 53[1] 54[1] 55[1] 56[1] 57[1] 58[1] 59[1] 60[1] 61[1] 62[1] 63[1] 64[1] 65[1] 66[1] 67[1] 68[1] 69[1] 70[1] 71[1] 72[1] 73[1] 74[1] 75[1] 76[1] 77[1] 78[1] 79[1] 80[1] 81[1] 82[1] 83[1] 84[1] 85[1] 86[1] 87# bind to tibble
hysteresis <- do.call(bind_rows, hysteresis_list) %>%
left_join(mysiteyrs %>% select(site_name, basin, subbasin, region, WaterYear, numdays, totalyield_site)) %>%
mutate(propyield = totalyieldevent_little / totalyield_site)Distribution of hystersis index by subbasin
hysteresis %>%
ggplot(aes(x = subbasin, y = hysteresis, fill = subbasin)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(height = 0, width = 0.25, alpha = 0.2) +
xlab("Sub-basin") + ylab("Xue hysteresis index") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.text.x = element_text(angle = 45, vjust = 0.5, hjust = 0.4))
Magnitude of hysteresis by length of event
hysteresis %>%
ggplot(aes(x = eventdays, y = abs(hysteresis), color = subbasin)) +
geom_point() +
xlab("Event length (days)") + ylab("Xue hysteresis index (absolute value)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
Magnitude of hysteresis by length of event, by subbasin and site
hysteresis %>%
#filter(eventdays >= 6) %>%
ggplot(aes(x = eventdays, y = abs(hysteresis), color = site_name, group = site_name)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap(~subbasin, scales = "free") +
xlab("Event length (days)") + ylab("Xue hysteresis index (absolute value)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
Magnitude of hysteresis by proportion of annual yield in event
hysteresis %>%
ggplot(aes(x = propyield, y = abs(hysteresis), color = subbasin)) +
geom_point() +
xlab("Proportion of annual yield in event") + ylab("Xue hysteresis index (absolute value)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
Magnitude of hysteresis by proportion of annual yield in event, by subbasin and site
hysteresis %>%
ggplot(aes(x = propyield, y = abs(hysteresis), color = site_name, group = site_name)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap(~subbasin, scales = "free") +
xlab("Proportion of annual yield in event") + ylab("Xue hysteresis index (absolute value)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
Summarize hysteresis by site and year
temp <- hysteresis %>%
group_by(subbasin, site_name, WaterYear) %>%
summarize(freq = n(),
meanpropyield = mean(propyield),
meanhyst = mean(abs(hysteresis))) %>%
mutate(logmeanhyst = log(meanhyst)) %>%
ungroup()
temp# A tibble: 87 × 7
subbasin site_name WaterYear freq meanpropyield meanhyst logmeanhyst
<chr> <chr> <dbl> <int> <dbl> <dbl> <dbl>
1 Big Creek Big Creek NWIS 2019 2 0.435 2.02 0.703
2 Big Creek Big Creek NWIS 2020 3 0.294 0.954 -0.0473
3 Big Creek Big Creek NWIS 2021 6 0.143 0.878 -0.130
4 Big Creek Big Creek NWIS 2022 7 0.134 0.746 -0.294
5 Big Creek Hallowat Creek… 2020 2 0.465 0.868 -0.141
6 Big Creek Hallowat Creek… 2022 6 0.162 1.05 0.0452
7 Big Creek NicolaCreek 2019 2 0.284 3.96 1.38
8 Big Creek WernerCreek 2019 2 0.407 1.84 0.608
9 Big Creek WernerCreek 2020 4 0.202 1.43 0.358
10 Coal Creek CycloneCreekUp… 2019 2 0.436 2.84 1.04
# ℹ 77 more rowsCharacterize sites/subbasins along axes of event frequency and magnitude (mean annual proportion of flow within event).
hull <- temp %>% group_by(subbasin) %>% slice(chull(freq, meanpropyield))
temp %>%
ggplot(aes(x = freq, y = meanpropyield)) +
#geom_polygon(data = hull, aes(fill = subbasin), alpha = 0.5) +
#geom_smooth(color = "black", method = "lm") +
geom_point(aes(color = subbasin)) +
xlab("Events per year") + ylab("Proportion of annual yield in event (mean over events)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
Plot the mean annual magnitude (absolute value) of event-specific hysteresis by event frequency (events per year), with LOESS smoothing.
p1 <- temp %>%
ggplot(aes(x = freq, y = (meanhyst))) +
geom_smooth(color = "black") +
#geom_smooth(color = "black", method = "glm", formula = y~x, method.args = list(family = gaussian(link = 'log'))) +
geom_point(aes(color = subbasin)) +
xlab("Events per year") + ylab("Mean annual magnitude of hysteresis") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p2 <- temp %>%
ggplot(aes(x = meanpropyield, y = (meanhyst))) +
geom_smooth(color = "black") +
#geom_smooth(color = "black", method = "glm", formula = y~x, method.args = list(family = gaussian(link = 'log'))) +
geom_point(aes(color = subbasin)) +
xlab("Proportion of annual yield in event") + ylab("Mean annual magnitude of hysteresis") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggarrange(p1, p2, nrow = 1, common.legend = TRUE, legend = "right")
Same as above but plot on log scale (magnitude):
p1 <- temp %>%
ggplot(aes(x = freq, y = logmeanhyst)) +
geom_smooth(color = "black", method = "lm") +
geom_point(aes(color = subbasin)) +
xlab("Events per year") + ylab("(log) Mean annual magnitude of hysteresis") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p2 <- temp %>%
ggplot(aes(x = meanpropyield, y = logmeanhyst)) +
geom_smooth(color = "black", method = "lm") +
geom_point(aes(color = subbasin)) +
xlab("Proportion of annual yield in event") + ylab("(log) Mean annual magnitude of hysteresis") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggarrange(p1, p2, nrow = 1, common.legend = TRUE, legend = "right")
Replace regression line with convex hulls
hull <- temp %>% group_by(subbasin) %>% slice(chull(freq, logmeanhyst))
p1 <- temp %>%
ggplot(aes(x = freq, y = logmeanhyst)) +
geom_polygon(data = hull, aes(fill = subbasin), alpha = 0.4) +
geom_point(aes(color = subbasin)) +
xlab("Events per year") + ylab("(log) Mean annual magnitude of hysteresis") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
hull <- temp %>% group_by(subbasin) %>% slice(chull(meanpropyield, logmeanhyst))
p2 <- temp %>%
ggplot(aes(x = meanpropyield, y = logmeanhyst)) +
geom_polygon(data = hull, aes(fill = subbasin), alpha = 0.4) +
geom_point(aes(color = subbasin)) +
xlab("Proportion of annual yield in event") + ylab("(log) Mean annual magnitude of hysteresis") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggarrange(p1, p2, nrow = 1, common.legend = TRUE, legend = "right")
Alternatively, plot frequency by magnitude and with point color/size. I don’t think this tells the story as well as the convex hull plots.
temp %>%
ggplot(aes(x = meanpropyield, y = freq)) +
# geom_polygon(data = hull, aes(fill = subbasin), alpha = 0.5) +
geom_point(aes(color = subbasin, size = logmeanhyst)) +
scale_size_continuous(range = c(0.1,5)) +
xlab("Proportion of annual yield in event") + ylab("Events per year") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
temp %>%
ggplot(aes(x = meanpropyield, y = freq)) +
# geom_polygon(data = hull, aes(fill = subbasin), alpha = 0.5) +
geom_point(aes(color = logmeanhyst), size = 2) +
xlab("Proportion of annual yield in event") + ylab("Events per year") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
Conduct baseflow separation, event delineation, and dynamic time warping as above, but only for a single site year. For troubleshooting and data visualization
Set parameters and site-year
alp <- 0.95
numpass <- 3
thresh <- 0.75
dd <- dat_clean_hyst %>% filter(site_name == "West Brook Reservoir", WaterYear == 2023)baseflow separation
dd <- dd %>%
filter(!is.na(Yield_mm_big)) %>%
mutate(bf = baseflowB(Yield_mm_big, alpha = alp, passes = numpass)$bf,
bfi = baseflowB(Yield_mm_big, alpha = alp, passes = numpass)$bfi) %>%
ungroup()
head(dd)Event delineation
events <- eventBaseflow(dd$Yield_mm_big, BFI_Th = thresh, bfi = dd$bfi)
events <- events %>% mutate(len = end - srt + 1)
(events)Tidy events: now add variables to the Big G time series data specifying events and non-events
# define positions of non-events
srt <- c(1)
end <- c(events$srt[1]-1)
for (i in 2:(dim(events)[1])) {
srt[i] <- events$end[i-1]+1
end[i] <- events$srt[i]-1
}
nonevents <- data.frame(tibble(srt, end) %>% mutate(len = end - srt) %>% filter(len >= 0) %>% select(-len) %>% add_row(srt = events$end[dim(events)[1]]+1, end = dim(dd)[1]))
# create vectors of binary event/non-event and event IDs
isevent_vec <- rep(2, times = dim(dd)[1])
eventid_vec <- rep(NA, times = dim(dd)[1])
for (i in 1:dim(events)[1]) {
isevent_vec[c(events[i,1]:events[i,2])] <- 1
eventid_vec[c(events[i,1]:events[i,2])] <- i
}
# create vector of non-event IDs
noneventid_vec <- rep(NA, times = dim(dd)[1])
for (i in 1:dim(nonevents)[1]) { noneventid_vec[c(nonevents[i,1]:nonevents[i,2])] <- i }
# create vector of "agnostic events": combined hydro events and non-events
agnevents <- rbind(events %>% select(srt, end) %>% mutate(event = 1), nonevents %>% mutate(event = 0)) %>% arrange((srt))
agneventid_vec <- c()
for (i in 1:dim(agnevents)[1]){ agneventid_vec[c(agnevents[i,1]:agnevents[i,2])] <- i }
# add event/non-event vectors to Big G data
dd <- dd %>%
mutate(isevent = isevent_vec,
eventid = eventid_vec,
noneventid = noneventid_vec,
agneventid = agneventid_vec,
big_event_yield = ifelse(isevent_vec == 1, Yield_mm_big, NA),
big_nonevent_yield = ifelse(isevent_vec == 2, Yield_mm_big, NA),
big_event_quick = big_event_yield - bf) %>%
rename(big_bf = bf, big_bfi = bfi)
(dd)View time series data with Big G event delineation
dd %>% select(date, Yield_mm_big, Yield_mm, big_event_yield, big_nonevent_yield) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "Yield (mm)") %>% dyOptions(fillGraph = TRUE, drawGrid = FALSE, axisLineWidth = 1.5)View hysteresis loops: g~G
dd %>%
filter(isevent == 1) %>%
ggplot(aes(x = logYield_big, y = logYield, color = date)) +
geom_segment(aes(xend = c(tail(logYield_big, n = -1), NA), yend = c(tail(logYield, n = -1), NA)),
arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_point() +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
scale_color_gradientn(colors = cet_pal(250, name = "r1"), trans = "date") +
facet_wrap(~eventid) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())Select and view specific event and normalize data
ddd <- dd %>% filter(eventid == 2)
dddd <- dd %>% filter(date %in% c(min(ddd$date)-1, max(ddd$date)+1))
ddd <- ddd %>%
bind_rows(dddd) %>%
arrange(date) %>%
mutate(weight = (logYield_big - min(logYield_big)) / (max(logYield_big) - min(logYield_big)),
yield_little_norm = (logYield - min(logYield)) / (max(logYield) - min(logYield)),
yield_big_norm = (logYield_big - min(logYield_big)) / (max(logYield_big) - min(logYield_big)))
ddd %>%
ggplot(aes(x = yield_big_norm, y = yield_little_norm, color = date)) +
geom_segment(aes(xend = c(tail(yield_big_norm, n = -1), NA), yend = c(tail(yield_little_norm, n = -1), NA)),
arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_point() +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
scale_color_gradientn(colors = cet_pal(250, name = "r1"), trans = "date") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())Align data using dynamic time warping
align <- dtw(x = unlist(ddd %>% select(yield_big_norm)), y = unlist(ddd %>% select(yield_little_norm)), step = asymmetric, keep = TRUE)
# align <- dtw(x = unlist(ddd %>% select(Yield_mm_big)), y = unlist(ddd %>% select(Yield_mm)), step = asymmetric, keep = TRUE)
plot(align, type = "threeway")
plot(align, type = "twoway", offset = -1)Calculate hysteresis index as in Xue et al (2024)
ddd <- ddd %>% mutate(distance = align$index1 - align$index2)
ddd %>%
ggplot(aes(x = date, y = distance*weight)) +
geom_bar(stat = "identity") +
theme_bw()
sum(ddd$weight * ddd$distance) / sum(ddd$weight)Suitability of existing modeling techniques
Purpose: Quantify the suitability of existing modeling techniques for predicting streamflow in headwater systems.
Approach:
Pull HUC-10 streamflow predictions from Montana Climate Office’s (University of Montana) Streamflow API
Get HUC-10 watershed codes
myhucs <- c()
for (i in 1:dim(siteinfo_sp)[1]) { myhucs[i] <- get_huc(AOI = siteinfo_sp[i,], type = "huc10")$huc10 }
siteinfo <- siteinfo %>% mutate(huc10 = myhucs)
myhucs <- unique(myhucs)
length(myhucs)[1] 16Query UM Climatology
request = httr::GET(
# can replace this with /predictions/raw, the only query parameter that isn't shared is aggregations.
"https://data.climate.umt.edu/streamflow-api/predictions/",
query = list(
locations = paste(myhucs, collapse = ","),
date_start = "2015-01-01",
date_end = "2025-01-01",
aggregations = "mean",
as_csv = TRUE,
units = "mm"
)
)
umpreds <- httr::content(request)
umpreds <- umpreds %>% rename(huc10 = location)
print(umpreds)# A tibble: 58,464 × 5
huc10 version date metric value
<chr> <chr> <date> <chr> <dbl>
1 0108020106 vPUB2025 2015-01-01 mean 2.14
2 0108020106 vPUB2025 2015-01-02 mean 1.86
3 0108020106 vPUB2025 2015-01-03 mean 1.93
4 0108020106 vPUB2025 2015-01-04 mean 4.00
5 0108020106 vPUB2025 2015-01-05 mean 4.04
6 0108020106 vPUB2025 2015-01-06 mean 3.07
7 0108020106 vPUB2025 2015-01-07 mean 2.43
8 0108020106 vPUB2025 2015-01-08 mean 2.16
9 0108020106 vPUB2025 2015-01-09 mean 2.03
10 0108020106 vPUB2025 2015-01-10 mean 1.90
# ℹ 58,454 more rowsView data
umpreds %>%
ggplot(aes(x = date, y = value)) +
geom_line() +
facet_wrap(~huc10, scales = "free_y")
Join to observed data
dat_umpred <- dat_clean %>%
left_join(siteinfo %>% select(site_name, huc10)) %>%
left_join(umpreds %>% select(huc10, date, value)) View example sites ::: panel-tabset #### Jimmy
dat_umpred %>% filter(site_name == "Jimmy Brook") %>% select(date, Yield_mm, value) %>% dygraph() %>% dySeries("Yield_mm", color = "black") %>% dySeries("value", color = "red") %>% dyRangeSelector() dat_umpred %>% filter(site_name == "Staunton River 06") %>% select(date, Yield_mm, value) %>% dygraph() %>% dySeries("Yield_mm", color = "black") %>% dySeries("value", color = "red") %>% dyRangeSelector() dat_umpred %>% filter(site_name == "Hallowat Creek NWIS") %>% select(date, Yield_mm, value) %>% dygraph() %>% dySeries("Yield_mm", color = "black") %>% dySeries("value", color = "red") %>% dyRangeSelector() dat_umpred %>% filter(site_name == "SF Spread Creek Lower") %>% select(date, Yield_mm, value) %>% dygraph() %>% dySeries("Yield_mm", color = "black") %>% dySeries("value", color = "red") %>% dyRangeSelector() dat_umpred %>% filter(site_name == "Dugout Creek") %>% select(date, Yield_mm, value) %>% dygraph() %>% dySeries("Yield_mm", color = "black") %>% dySeries("value", color = "red") %>% dyRangeSelector() dat_umpred %>% filter(site_name == "Donner Blitzen ab Indian NWIS") %>% select(date, Yield_mm, value) %>% dygraph() %>% dySeries("Yield_mm", color = "black") %>% dySeries("value", color = "red") %>% dyRangeSelector() :::
NSE ratings per Moriasi et al. “Model evaluation guidelines for systematic quantification of accuracy in watershed simulations.” Transactions of the ASABE 50.3 (2007): 885-900.
Note that there is no effort to ensure consistant data availability among sites/basins/years. Would it be better to restrict these?
# all data...
dat_umpred %>%
group_by(site_name, subbasin) %>%
summarize(nse = NSE(sim = log(value), obs = log(Yield_mm))) %>%
ggplot(aes(x = subbasin, y = nse)) +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0.50, alpha = 0.3, fill = "red") +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.50, ymax = 0.65, alpha = 0.3, fill = "orange") +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.65, ymax = 0.75, alpha = 0.3, fill = "yellow") +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.75, ymax = 1.00, alpha = 0.3, fill = "green") +
geom_boxplot(fill = "grey", outlier.shape = NA) +
geom_jitter(height = 0, width = 0.2, shape = 1) +
geom_abline(slope = 0, intercept = 1, color = "black", linetype = "dashed") +
xlab("Sub-basin") + ylab("Nash-Sutcliffe efficiency") + ggtitle("All data") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.text.x = element_text(angle = 45, vjust = 0.5, hjust = 0.4)) +
scale_y_continuous(limits = c(-8.5,1), expand = c(0,0))
# summer only
dat_umpred %>%
filter(Month %in% c(7:9)) %>%
group_by(site_name, subbasin) %>%
summarize(nse = NSE(sim = log(value), obs = log(Yield_mm))) %>%
ggplot(aes(x = subbasin, y = nse)) +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0.50, alpha = 0.3, fill = "red") +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.50, ymax = 0.65, alpha = 0.3, fill = "orange") +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.65, ymax = 0.75, alpha = 0.3, fill = "yellow") +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.75, ymax = 1.00, alpha = 0.3, fill = "green") +
geom_boxplot(fill = "grey", outlier.shape = NA) +
geom_jitter(height = 0, width = 0.2, shape = 1) +
geom_abline(slope = 0, intercept = 1, color = "black", linetype = "dashed") +
xlab("Sub-basin") + ylab("Nash-Sutcliffe efficiency") + ggtitle("July-September") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.text.x = element_text(angle = 45, vjust = 0.5, hjust = 0.4)) +
scale_y_continuous(limits = c(-12.5,1), expand = c(0,0))
How does NSE vary through time?
# all data...
dat_umpred %>%
group_by(site_name, subbasin, Month) %>%
summarize(nse = NSE(sim = log(value), obs = log(Yield_mm))) %>%
ggplot(aes(x = Month, y = nse)) +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0.50, alpha = 0.3, fill = "red") +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.50, ymax = 0.65, alpha = 0.3, fill = "orange") +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.65, ymax = 0.75, alpha = 0.3, fill = "yellow") +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.75, ymax = 1.00, alpha = 0.3, fill = "green") +
geom_point(aes(group = site_name)) +
geom_abline(slope = 0, intercept = 1, color = "black", linetype = "dashed") +
xlab("Month") + ylab("Nash-Sutcliffe efficiency") + ggtitle("All data") +
facet_wrap(~subbasin, scales = "free_y") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
scale_y_continuous(expand = c(0,0))
# West Brook
dat_umpred %>%
mutate(yearmonth = floor_date(date, "month")) %>%
filter(subbasin == "West Brook") %>%
group_by(site_name, subbasin, yearmonth) %>%
summarize(nobs = n(), nse = NSE(sim = log(value), obs = log(Yield_mm))) %>%
ungroup() %>%
filter(nobs >= 25) %>%
ggplot(aes(x = yearmonth, y = nse)) +
# annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0.50, alpha = 0.3, fill = "red") +
# annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.50, ymax = 0.65, alpha = 0.3, fill = "orange") +
# annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.65, ymax = 0.75, alpha = 0.3, fill = "yellow") +
# annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.75, ymax = 1.00, alpha = 0.3, fill = "green") +
geom_point(aes(color = site_name)) +
geom_line(aes(color = site_name)) +
geom_abline(slope = 0, intercept = 1, color = "black", linetype = "dashed") +
xlab("Time") + ylab("Nash-Sutcliffe efficiency") + ggtitle("West Brook (truncated y-axis limits)") +
#facet_wrap(~subbasin, scales = "free") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
ylim(-50,1)
daily <- dat_umpred %>%
group_by(site_name, subbasin) %>%
summarize(nse = NSE(sim = log(value), obs = log(Yield_mm))) %>%
ungroup() %>%
mutate(timescale = "daily")
weekly <- dat_umpred %>%
mutate(date = floor_date(date, "week")) %>%
group_by(site_name, subbasin, date) %>%
summarise(nobs = n(), Yield_mm = mean(Yield_mm), value = mean(value)) %>%
ungroup() %>%
filter(nobs == 7) %>%
group_by(site_name, subbasin) %>%
summarize(nse = NSE(sim = log(value), obs = log(Yield_mm))) %>%
ungroup() %>%
mutate(timescale = "weekly")
monthly <- dat_umpred %>%
mutate(date = floor_date(date, "month")) %>%
group_by(site_name, subbasin, date) %>%
summarise(nobs = n(), Yield_mm = mean(Yield_mm), value = mean(value)) %>%
ungroup() %>%
filter(nobs >= 27) %>%
group_by(site_name, subbasin) %>%
summarize(nse = NSE(sim = log(value), obs = log(Yield_mm))) %>%
ungroup() %>%
mutate(timescale = "monthly")
timescale <- bind_rows(daily, weekly, monthly) %>%
mutate(timescale = factor(timescale, levels = c("daily", "weekly", "monthly")))
timescale %>%
ggplot(aes(x = timescale, y = nse, group = site_name)) +
geom_point() +
geom_line() +
xlab("Time scale") + ylab("Nash-Sutcliffe efficiency") +
theme_bw() + ylim(-6,1) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
timescale %>%
ggplot(aes(x = timescale, y = nse, group = site_name)) +
geom_point() +
geom_line() +
xlab("Time scale") + ylab("Nash-Sutcliffe efficiency") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
facet_wrap(~subbasin, scales = "free_y")
Purpose: Demonstrate spatial variation in drought-related low flow conditions across headwater stream networks.
Approach:
Generate fixed and time varying (by day of year) drought/low flow thresholds from long-term (1970-2025) big G streamflow data. Quantiles in ~0.05 increments from 0.02 to 0.50.
# calculate site-level fixed low flow thresholds
dat_clean_big_fixed <- dat_clean_big %>%
filter(!is.na(basin)) %>%
group_by(site_name, basin, subbasin, region) %>%
summarize(thresh_50_fix = quantile(Yield_mm, probs = 0.50, na.rm = TRUE),
thresh_45_fix = quantile(Yield_mm, probs = 0.45, na.rm = TRUE),
thresh_40_fix = quantile(Yield_mm, probs = 0.40, na.rm = TRUE),
thresh_35_fix = quantile(Yield_mm, probs = 0.35, na.rm = TRUE),
thresh_30_fix = quantile(Yield_mm, probs = 0.30, na.rm = TRUE),
thresh_25_fix = quantile(Yield_mm, probs = 0.25, na.rm = TRUE),
thresh_20_fix = quantile(Yield_mm, probs = 0.20, na.rm = TRUE),
thresh_15_fix = quantile(Yield_mm, probs = 0.15, na.rm = TRUE),
thresh_10_fix = quantile(Yield_mm, probs = 0.10, na.rm = TRUE),
thresh_05_fix = quantile(Yield_mm, probs = 0.05, na.rm = TRUE),
thresh_02_fix = quantile(Yield_mm, probs = 0.02, na.rm = TRUE)) %>%
ungroup()
(dat_clean_big_fixed)# A tibble: 8 × 15
site_name basin subbasin region thresh_50_fix thresh_45_fix thresh_40_fix
<chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
1 Battle Run NW… Pine… Piney R… Shen 0.586 0.513 0.439
2 Donner Blitze… Donn… Donner … Oreg 0.272 0.246 0.227
3 North Fork Fl… Flat… Flathead Flat 0.729 0.644 0.579
4 Pacific Creek… Snak… Snake R… Snake 0.364 0.336 0.313
5 Rapidan River… Stau… Staunto… Shen 0.830 0.731 0.642
6 South River C… West… West Br… Mass 1.34 1.18 1.02
7 South River H… Pain… Paine R… Shen 0.704 0.624 0.557
8 Yellowstone R… Shie… Shields… Shiel… 0.500 0.466 0.438
# ℹ 8 more variables: thresh_35_fix <dbl>, thresh_30_fix <dbl>,
# thresh_25_fix <dbl>, thresh_20_fix <dbl>, thresh_15_fix <dbl>,
# thresh_10_fix <dbl>, thresh_05_fix <dbl>, thresh_02_fix <dbl># calculate site-level variable (by doy) low flow thresholds
dat_clean_big_variable <- dat_clean_big %>%
filter(!is.na(basin)) %>%
group_by(site_name, basin, subbasin, region, doy_calendar) %>%
summarize(thresh_50_var = quantile(Yield_mm, probs = 0.50, na.rm = TRUE),
thresh_45_var = quantile(Yield_mm, probs = 0.45, na.rm = TRUE),
thresh_40_var = quantile(Yield_mm, probs = 0.40, na.rm = TRUE),
thresh_35_var = quantile(Yield_mm, probs = 0.35, na.rm = TRUE),
thresh_30_var = quantile(Yield_mm, probs = 0.30, na.rm = TRUE),
thresh_25_var = quantile(Yield_mm, probs = 0.25, na.rm = TRUE),
thresh_20_var = quantile(Yield_mm, probs = 0.20, na.rm = TRUE),
thresh_15_var = quantile(Yield_mm, probs = 0.15, na.rm = TRUE),
thresh_10_var = quantile(Yield_mm, probs = 0.10, na.rm = TRUE),
thresh_05_var = quantile(Yield_mm, probs = 0.05, na.rm = TRUE),
thresh_02_var = quantile(Yield_mm, probs = 0.02, na.rm = TRUE)) %>%
ungroup()
(dat_clean_big_variable)# A tibble: 2,928 × 16
site_name basin subbasin region doy_calendar thresh_50_var thresh_45_var
<chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
1 Battle Run NW… Pine… Piney R… Shen 1 0.732 0.687
2 Battle Run NW… Pine… Piney R… Shen 2 0.754 0.702
3 Battle Run NW… Pine… Piney R… Shen 3 0.815 0.727
4 Battle Run NW… Pine… Piney R… Shen 4 0.860 0.802
5 Battle Run NW… Pine… Piney R… Shen 5 0.805 0.732
6 Battle Run NW… Pine… Piney R… Shen 6 0.751 0.690
7 Battle Run NW… Pine… Piney R… Shen 7 0.859 0.790
8 Battle Run NW… Pine… Piney R… Shen 8 0.787 0.732
9 Battle Run NW… Pine… Piney R… Shen 9 0.805 0.784
10 Battle Run NW… Pine… Piney R… Shen 10 0.820 0.732
# ℹ 2,918 more rows
# ℹ 9 more variables: thresh_40_var <dbl>, thresh_35_var <dbl>,
# thresh_30_var <dbl>, thresh_25_var <dbl>, thresh_20_var <dbl>,
# thresh_15_var <dbl>, thresh_10_var <dbl>, thresh_05_var <dbl>,
# thresh_02_var <dbl>Join drought thresholds derived from big G to little g’s
mydroughtlevels <- c("none", "q50", "q45", "q40", "q35", "q30", "q25", "q20", "q15", "q10", "q05", "q02")
# little gs
dat_clean_little_join <- dat_clean %>%
left_join(dat_clean_big_fixed %>% select(-c(site_name, subbasin))) %>%
left_join(dat_clean_big_variable %>% select(-c(site_name, subbasin))) %>%
mutate(month = month(date),
year = year(date),
designation = "little",
drought_fix = ifelse(Yield_mm <= thresh_50_fix & Yield_mm > thresh_45_fix, "q50",
ifelse(Yield_mm <= thresh_45_fix & Yield_mm > thresh_40_fix, "q45",
ifelse(Yield_mm <= thresh_40_fix & Yield_mm > thresh_35_fix, "q40",
ifelse(Yield_mm <= thresh_35_fix & Yield_mm > thresh_30_fix, "q35",
ifelse(Yield_mm <= thresh_30_fix & Yield_mm > thresh_25_fix, "q30",
ifelse(Yield_mm <= thresh_25_fix & Yield_mm > thresh_20_fix, "q25",
ifelse(Yield_mm <= thresh_20_fix & Yield_mm > thresh_15_fix, "q20",
ifelse(Yield_mm <= thresh_15_fix & Yield_mm > thresh_10_fix, "q15",
ifelse(Yield_mm <= thresh_10_fix & Yield_mm > thresh_05_fix, "q10",
ifelse(Yield_mm <= thresh_05_fix & Yield_mm > thresh_02_fix, "q05",
ifelse(Yield_mm <= thresh_02_fix, "q02", "none"))))))))))),
drought_var = ifelse(Yield_mm <= thresh_50_var & Yield_mm > thresh_45_var, "q50",
ifelse(Yield_mm <= thresh_45_var & Yield_mm > thresh_40_var, "q45",
ifelse(Yield_mm <= thresh_40_var & Yield_mm > thresh_35_var, "q40",
ifelse(Yield_mm <= thresh_35_var & Yield_mm > thresh_30_var, "q35",
ifelse(Yield_mm <= thresh_30_var & Yield_mm > thresh_25_var, "q30",
ifelse(Yield_mm <= thresh_25_var & Yield_mm > thresh_20_var, "q25",
ifelse(Yield_mm <= thresh_20_var & Yield_mm > thresh_15_var, "q20",
ifelse(Yield_mm <= thresh_15_var & Yield_mm > thresh_10_var, "q15",
ifelse(Yield_mm <= thresh_10_var & Yield_mm > thresh_05_var, "q10",
ifelse(Yield_mm <= thresh_05_var & Yield_mm > thresh_02_var, "q05",
ifelse(Yield_mm <= thresh_02_var, "q02", "none")))))))))))) %>%
mutate(drought_fix = factor(ifelse(is.na(Yield_mm), NA, drought_fix), levels = mydroughtlevels),
drought_var = factor(ifelse(is.na(Yield_mm), NA, drought_var), levels = mydroughtlevels))
# big gs
dat_clean_big_join <- dat_clean_big %>%
left_join(dat_clean_big_fixed %>% select(-c(site_name, subbasin))) %>%
left_join(dat_clean_big_variable %>% select(-c(site_name, subbasin))) %>%
mutate(month = month(date),
year = year(date),
designation = "big",
drought_fix = ifelse(Yield_mm <= thresh_50_fix & Yield_mm > thresh_45_fix, "q50",
ifelse(Yield_mm <= thresh_45_fix & Yield_mm > thresh_40_fix, "q45",
ifelse(Yield_mm <= thresh_40_fix & Yield_mm > thresh_35_fix, "q40",
ifelse(Yield_mm <= thresh_35_fix & Yield_mm > thresh_30_fix, "q35",
ifelse(Yield_mm <= thresh_30_fix & Yield_mm > thresh_25_fix, "q30",
ifelse(Yield_mm <= thresh_25_fix & Yield_mm > thresh_20_fix, "q25",
ifelse(Yield_mm <= thresh_20_fix & Yield_mm > thresh_15_fix, "q20",
ifelse(Yield_mm <= thresh_15_fix & Yield_mm > thresh_10_fix, "q15",
ifelse(Yield_mm <= thresh_10_fix & Yield_mm > thresh_05_fix, "q10",
ifelse(Yield_mm <= thresh_05_fix & Yield_mm > thresh_02_fix, "q05",
ifelse(Yield_mm <= thresh_02_fix, "q02", "none"))))))))))),
drought_var = ifelse(Yield_mm <= thresh_50_var & Yield_mm > thresh_45_var, "q50",
ifelse(Yield_mm <= thresh_45_var & Yield_mm > thresh_40_var, "q45",
ifelse(Yield_mm <= thresh_40_var & Yield_mm > thresh_35_var, "q40",
ifelse(Yield_mm <= thresh_35_var & Yield_mm > thresh_30_var, "q35",
ifelse(Yield_mm <= thresh_30_var & Yield_mm > thresh_25_var, "q30",
ifelse(Yield_mm <= thresh_25_var & Yield_mm > thresh_20_var, "q25",
ifelse(Yield_mm <= thresh_20_var & Yield_mm > thresh_15_var, "q20",
ifelse(Yield_mm <= thresh_15_var & Yield_mm > thresh_10_var, "q15",
ifelse(Yield_mm <= thresh_10_var & Yield_mm > thresh_05_var, "q10",
ifelse(Yield_mm <= thresh_05_var & Yield_mm > thresh_02_var, "q05",
ifelse(Yield_mm <= thresh_02_var, "q02", "none")))))))))))) %>%
mutate(drought_fix = factor(ifelse(is.na(Yield_mm), NA, drought_fix), levels = mydroughtlevels),
drought_var = factor(ifelse(is.na(Yield_mm), NA, drought_var), levels = mydroughtlevels))
dat_clean_big_join_sub <- dat_clean_big_join %>% left_join(yrdf) %>% filter(ecodyr == "yes")
# join data
dat_clean_join <- bind_rows(dat_clean_little_join, dat_clean_big_join_sub)Drought/low flow delineation is somewhat complicated by the fact that some streams simply have greater yield than others. For example, at groundwater-dominated sites, the above approach will never detect low flow conditions based on big G thresholds, but this doesn’t mean that flow at that site isn’t lower than normal (for that site). This is most obvious in the Snake River basin, where NF Spread Creek Upper never experiences drought (because this is presumably a gaining reach) and Rock and Grouse Creeks are in a perpetual state of drought (presumable these are losing reaches). This is a classic “At which level of organization do I standardize my data?” question: are general differences in flow volume among sites signal or noise? But perhaps more importantly, this is a question of “what is drought?” Is drought relative to some larger regional metric (e.g., big G)? Or is it a local phenomenon, where the specifics of individual streams and reaches matter.
For each site individually, generate (fixed) drought/low flow thresholds using the same quantiles same as above: ~0.05 increments from 0.02 to 0.50. Restrict data to selected basins, sites, and years with (nearly) complete summer (July, August, September) data over the selected periods/locations. (Standardization needs to be done over comparable time periods, at least among sites within basins).
Require 95% data availability across all water years for site to be included!
Organize data, get site-level low flow threshold values, and denote drought periods
# Require 95% data availability!
monthss <- c(7:9)
# grab data and bind, z-score Yield
dat_clean_sub <- bind_rows(
dat_clean %>% filter(subbasin == "West Brook", WaterYear %in% c(2020:2022), Month %in% monthss, !site_name %in% c("Mitchell Brook")) %>%
bind_rows(dat_clean_big %>% filter(subbasin == "West Brook", WaterYear %in% c(2020:2022), Month %in% monthss)),
dat_clean %>% filter(subbasin == "Staunton River", WaterYear %in% c(2019:2022), Month %in% monthss) %>%
bind_rows(dat_clean_big %>% filter(subbasin == "Staunton River", WaterYear %in% c(2019:2022), Month %in% monthss)),
dat_clean %>% filter(subbasin == "Big Creek", WaterYear %in% c(2019:2021), Month %in% monthss, !site_name %in% c("BigCreekLower", "LangfordCreekUpper", "NicolaCreek", "WernerCreek")) %>%
bind_rows(dat_clean_big %>% filter(basin == "Flathead", WaterYear %in% c(2019:2021), Month %in% monthss) %>% mutate(subbasin = "Big Creek")),
dat_clean %>% filter(subbasin == "Coal Creek", WaterYear %in% c(2019:2021), Month %in% monthss, !site_name %in% c("CycloneCreekMiddle", "CoalCreekMiddle", "CoalCreekHeadwaters")) %>%
bind_rows(dat_clean_big %>% filter(basin == "Flathead", WaterYear %in% c(2019:2021), Month %in% monthss) %>% mutate(subbasin = "Coal Creek")),
dat_clean %>% filter(subbasin == "Snake River", WaterYear %in% c(2020:2022), Month %in% monthss, !site_name %in% c("Spread Creek Dam")) %>%
bind_rows(dat_clean_big %>% filter(subbasin == "Snake River", WaterYear %in% c(2020:2022), Month %in% monthss)),
dat_clean %>% filter(subbasin == "Shields River", WaterYear %in% c(2019, 2020, 2023), Month %in% monthss, !site_name %in% c("Shields River Valley Ranch")) %>% group_by(site_name) %>%
bind_rows(dat_clean_big %>% filter(basin == "Shields River", WaterYear %in% c(2019, 2020, 2023), Month %in% monthss)),
dat_clean %>% filter(subbasin == "Duck Creek", WaterYear %in% c(2017:2022), Month %in% monthss) %>%
bind_rows(dat_clean_big %>% filter(basin == "Shields River", WaterYear %in% c(2017:2022), Month %in% monthss) %>% mutate(subbasin = "Duck Creek")),
dat_clean %>% filter(subbasin == "Donner Blitzen", WaterYear %in% c(2020:2022), Month %in% monthss) %>%
bind_rows(dat_clean_big %>% filter(basin == "Donner Blitzen", WaterYear %in% c(2020:2022), Month %in% monthss))
) %>%
group_by(site_name) %>%
mutate(z_Yield_mm = scale(Yield_mm, center = TRUE, scale = TRUE)[,1]) %>%
ungroup()
# get low flow thresholds
dat_clean_sub_thresh <- dat_clean_sub %>%
group_by(site_name) %>%
summarize(thresh_50_fix = quantile(z_Yield_mm, probs = 0.50, na.rm = TRUE),
thresh_45_fix = quantile(z_Yield_mm, probs = 0.45, na.rm = TRUE),
thresh_40_fix = quantile(z_Yield_mm, probs = 0.40, na.rm = TRUE),
thresh_35_fix = quantile(z_Yield_mm, probs = 0.35, na.rm = TRUE),
thresh_30_fix = quantile(z_Yield_mm, probs = 0.30, na.rm = TRUE),
thresh_25_fix = quantile(z_Yield_mm, probs = 0.25, na.rm = TRUE),
thresh_20_fix = quantile(z_Yield_mm, probs = 0.20, na.rm = TRUE),
thresh_15_fix = quantile(z_Yield_mm, probs = 0.15, na.rm = TRUE),
thresh_10_fix = quantile(z_Yield_mm, probs = 0.10, na.rm = TRUE),
thresh_05_fix = quantile(z_Yield_mm, probs = 0.05, na.rm = TRUE),
thresh_02_fix = quantile(z_Yield_mm, probs = 0.02, na.rm = TRUE)) %>%
ungroup()
dat_clean_sub_thresh# A tibble: 51 × 12
site_name thresh_50_fix thresh_45_fix thresh_40_fix thresh_35_fix
<chr> <dbl> <dbl> <dbl> <dbl>
1 Avery Brook -0.346 -0.361 -0.373 -0.380
2 Big Creek NWIS -0.380 -0.421 -0.474 -0.532
3 BigCreekUpper -0.435 -0.488 -0.532 -0.576
4 Buck Creek -0.339 -0.350 -0.369 -0.423
5 CoalCreekLower -0.378 -0.443 -0.507 -0.553
6 CoalCreekNorth -0.349 -0.429 -0.488 -0.548
7 Crandall Creek -0.319 -0.361 -0.390 -0.409
8 CycloneCreekLower -0.225 -0.313 -0.405 -0.457
9 CycloneCreekUpper -0.504 -0.555 -0.596 -0.617
10 Deep Creek -0.346 -0.380 -0.406 -0.442
# ℹ 41 more rows
# ℹ 7 more variables: thresh_30_fix <dbl>, thresh_25_fix <dbl>,
# thresh_20_fix <dbl>, thresh_15_fix <dbl>, thresh_10_fix <dbl>,
# thresh_05_fix <dbl>, thresh_02_fix <dbl># join thresholds to data and denote drought periods
dat_clean_sub <- dat_clean_sub %>%
left_join(dat_clean_sub_thresh) %>%
mutate(month = month(date),
year = year(date),
drought_fix = ifelse(z_Yield_mm <= thresh_50_fix & z_Yield_mm > thresh_45_fix, "q50",
ifelse(z_Yield_mm <= thresh_45_fix & z_Yield_mm > thresh_40_fix, "q45",
ifelse(z_Yield_mm <= thresh_40_fix & z_Yield_mm > thresh_35_fix, "q40",
ifelse(z_Yield_mm <= thresh_35_fix & z_Yield_mm > thresh_30_fix, "q35",
ifelse(z_Yield_mm <= thresh_30_fix & z_Yield_mm > thresh_25_fix, "q30",
ifelse(z_Yield_mm <= thresh_25_fix & z_Yield_mm > thresh_20_fix, "q25",
ifelse(z_Yield_mm <= thresh_20_fix & z_Yield_mm > thresh_15_fix, "q20",
ifelse(z_Yield_mm <= thresh_15_fix & z_Yield_mm > thresh_10_fix, "q15",
ifelse(z_Yield_mm <= thresh_10_fix & z_Yield_mm > thresh_05_fix, "q10",
ifelse(z_Yield_mm <= thresh_05_fix & z_Yield_mm > thresh_02_fix, "q05",
ifelse(z_Yield_mm <= thresh_02_fix, "q02", "none")))))))))))) %>%
mutate(drought_fix = factor(ifelse(is.na(Yield_mm), NA, drought_fix), levels = mydroughtlevels))
dat_clean_sub# A tibble: 16,243 × 31
site_name basin subbasin region date flow_mean tempc_mean Yield_mm
<chr> <chr> <chr> <chr> <date> <dbl> <dbl> <dbl>
1 Avery Brook West Br… West Br… Mass 2020-07-01 1.34 16.0 0.446
2 Avery Brook West Br… West Br… Mass 2020-07-02 0.963 16.1 0.321
3 Avery Brook West Br… West Br… Mass 2020-07-03 2.40 17.3 0.800
4 Avery Brook West Br… West Br… Mass 2020-07-04 3.31 17.4 1.10
5 Avery Brook West Br… West Br… Mass 2020-07-05 1.38 17.5 0.460
6 Avery Brook West Br… West Br… Mass 2020-07-06 0.965 17.5 0.322
7 Avery Brook West Br… West Br… Mass 2020-07-07 0.778 17.3 0.259
8 Avery Brook West Br… West Br… Mass 2020-07-08 0.795 17.3 0.265
9 Avery Brook West Br… West Br… Mass 2020-07-09 1.07 17.9 0.356
10 Avery Brook West Br… West Br… Mass 2020-07-10 14.3 19.4 4.76
# ℹ 16,233 more rows
# ℹ 23 more variables: CalendarYear <dbl>, Month <dbl>, MonthName <fct>,
# WaterYear <dbl>, DayofYear <dbl>, logYield <dbl>, designation <chr>,
# doy_calendar <dbl>, z_Yield_mm <dbl>, thresh_50_fix <dbl>,
# thresh_45_fix <dbl>, thresh_40_fix <dbl>, thresh_35_fix <dbl>,
# thresh_30_fix <dbl>, thresh_25_fix <dbl>, thresh_20_fix <dbl>,
# thresh_15_fix <dbl>, thresh_10_fix <dbl>, thresh_05_fix <dbl>, …Create heatmap functions
# fixed drought threshold
heatmapfun_fix <- function(subbas, months, bigG) {
dd <- dat_clean_join %>% filter(subbasin == subbas | site_name == bigG, month %in% months)
mysites <- c(unique(unlist(dd %>% filter(designation == "big") %>% select(site_name))),
unique(unlist(dd %>% filter(designation == "little") %>% select(site_name))))
myrect <- dd %>% group_by(WaterYear) %>% summarize(mindate = min(date), maxdate = max(date)) %>% ungroup()
p <- dd %>%
ggplot() +
geom_tile(aes(x = date, y = factor(site_name, levels = rev(mysites)), fill = drought_fix)) +
scale_fill_viridis(option = "A", direction = -1, discrete = TRUE, limits = mydroughtlevels) +
geom_rect(data = myrect, aes(xmin = mindate, xmax = maxdate, ymin = length(mysites)-0.5, ymax = length(mysites)+0.5),
color = "grey70", fill = NA, size = 1.25) +
xlab("Date") + ylab("Site") +
#facet_wrap(~WaterYear, scales = "free_x") +
facet_wrap2(~WaterYear, scales = "free_x", nrow = 3, ncol = 3, trim_blank = FALSE) +
theme_bw() + theme(axis.title = element_blank())
return(p)
}
# variable drought threshold
heatmapfun_var <- function(subbas, months, bigG) {
dd <- dat_clean_join %>% filter(subbasin == subbas | site_name == bigG, month %in% months)
mysites <- c(unique(unlist(dd %>% filter(designation == "big") %>% select(site_name))),
unique(unlist(dd %>% filter(designation == "little") %>% select(site_name))))
myrect <- dd %>% group_by(WaterYear) %>% summarize(mindate = min(date), maxdate = max(date)) %>% ungroup()
p <- dd %>%
ggplot() +
geom_tile(aes(x = date, y = factor(site_name, levels = rev(mysites)), fill = drought_var)) +
scale_fill_viridis(option = "A", direction = -1, discrete = TRUE, limits = mydroughtlevels) +
geom_rect(data = myrect, aes(xmin = mindate, xmax = maxdate, ymin = length(mysites)-0.5, ymax = length(mysites)+0.5),
color = "grey70", fill = NA, size = 1.25) +
xlab("Date") + ylab("Site") +
#facet_wrap(~WaterYear, scales = "free_x") +
facet_wrap2(~WaterYear, scales = "free_x", nrow = 3, ncol = 3, trim_blank = FALSE) +
theme_bw() + theme(axis.title = element_blank())
return(p)
}
# site-level drought threshold
heatmapfun_site <- function(subbas, months, bigG) {
dd <- dat_clean_sub %>% filter(subbasin == subbas)
mysites <- c(unique(unlist(dd %>% filter(site_name != bigG) %>% select(site_name))), bigG)
myrect <- dd %>% group_by(WaterYear) %>% summarize(mindate = min(date), maxdate = max(date)) %>% ungroup()
p <- dd %>%
ggplot() +
geom_tile(aes(x = date, y = factor(site_name, levels = (mysites)), fill = drought_fix)) +
scale_fill_viridis(option = "A", direction = -1, discrete = TRUE, limits = mydroughtlevels) +
geom_rect(data = myrect, aes(xmin = mindate, xmax = maxdate, ymin = length(mysites)-0.5, ymax = length(mysites)+0.5),
color = "grey70", fill = NA, size = 1.25) +
xlab("Date") + ylab("Site") +
#facet_wrap(~WaterYear, scales = "free_x") +
facet_wrap2(~WaterYear, scales = "free_x", nrow = 2, ncol = 3, trim_blank = FALSE) +
theme_bw() + theme(axis.title = element_blank())
return(p)
}heatmapfun_fix(subbas = "West Brook", bigG = "South River Conway NWIS", months = c(7:9))
heatmapfun_fix(subbas = "Paine Run", bigG = "South River Harriston NWIS", months = c(7:9))
heatmapfun_fix(subbas = "Staunton River", bigG = "Rapidan River NWIS", months = c(7:9))
heatmapfun_fix(subbas = "Big Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9))
heatmapfun_fix(subbas = "Coal Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9))
heatmapfun_fix(subbas = "McGee Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9))
heatmapfun_fix(subbas = "Snake River", bigG = "Pacific Creek at Moran NWIS", months = c(7:9))
heatmapfun_fix(subbas = "Shields River", bigG = "Yellowstone River Livingston NWIS", months = c(7:9))
heatmapfun_fix(subbas = "Duck Creek", bigG = "Yellowstone River Livingston NWIS", months = c(7:9))
heatmapfun_fix(subbas = "Donner Blitzen", bigG = "Donner Blitzen River nr Frenchglen NWIS", months = c(7:9))
heatmapfun_var(subbas = "West Brook", bigG = "South River Conway NWIS", months = c(7:9))
heatmapfun_var(subbas = "Paine Run", bigG = "South River Harriston NWIS", months = c(7:9))
heatmapfun_var(subbas = "Staunton River", bigG = "Rapidan River NWIS", months = c(7:9))
heatmapfun_var(subbas = "Big Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9))
heatmapfun_var(subbas = "Coal Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9))
heatmapfun_var(subbas = "McGee Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9))
heatmapfun_var(subbas = "Snake River", bigG = "Pacific Creek at Moran NWIS", months = c(7:9))
heatmapfun_var(subbas = "Shields River", bigG = "Yellowstone River Livingston NWIS", months = c(7:9))
heatmapfun_var(subbas = "Duck Creek", bigG = "Yellowstone River Livingston NWIS", months = c(7:9))
heatmapfun_var(subbas = "Donner Blitzen", bigG = "Donner Blitzen River nr Frenchglen NWIS", months = c(7:9))
heatmapfun_site(subbas = "West Brook", bigG = "South River Conway NWIS", months = c(7:9))
heatmapfun_site(subbas = "Staunton River", bigG = "Rapidan River NWIS", months = c(7:9))
heatmapfun_site(subbas = "Big Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9))
heatmapfun_site(subbas = "Coal Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9))
heatmapfun_site(subbas = "Snake River", bigG = "Pacific Creek at Moran NWIS", months = c(7:9))
heatmapfun_site(subbas = "Shields River", bigG = "Yellowstone River Livingston NWIS", months = c(7:9))
heatmapfun_site(subbas = "Duck Creek", bigG = "Yellowstone River Livingston NWIS", months = c(7:9))
heatmapfun_site(subbas = "Donner Blitzen", bigG = "Donner Blitzen River nr Frenchglen NWIS", months = c(7:9))
Calculate daily deficit, then summarize deficit magnitude and duration by site and month.
# calculate daily deficit
dat_clean_join_deficit <- dat_clean_join %>%
mutate(deficit_50_fix = ifelse(Yield_mm < thresh_50_fix, thresh_50_fix - Yield_mm, 0),
deficit_45_fix = ifelse(Yield_mm < thresh_45_fix, thresh_45_fix - Yield_mm, 0),
deficit_40_fix = ifelse(Yield_mm < thresh_40_fix, thresh_40_fix - Yield_mm, 0),
deficit_35_fix = ifelse(Yield_mm < thresh_35_fix, thresh_35_fix - Yield_mm, 0),
deficit_30_fix = ifelse(Yield_mm < thresh_30_fix, thresh_30_fix - Yield_mm, 0),
deficit_25_fix = ifelse(Yield_mm < thresh_25_fix, thresh_25_fix - Yield_mm, 0),
deficit_20_fix = ifelse(Yield_mm < thresh_20_fix, thresh_20_fix - Yield_mm, 0),
deficit_15_fix = ifelse(Yield_mm < thresh_15_fix, thresh_15_fix - Yield_mm, 0),
deficit_10_fix = ifelse(Yield_mm < thresh_10_fix, thresh_10_fix - Yield_mm, 0),
deficit_05_fix = ifelse(Yield_mm < thresh_05_fix, thresh_05_fix - Yield_mm, 0),
deficit_02_fix = ifelse(Yield_mm < thresh_05_fix, thresh_05_fix - Yield_mm, 0),
deficit_50_var = ifelse(Yield_mm < thresh_50_var, thresh_50_var - Yield_mm, 0),
deficit_45_var = ifelse(Yield_mm < thresh_45_var, thresh_45_var - Yield_mm, 0),
deficit_40_var = ifelse(Yield_mm < thresh_40_var, thresh_40_var - Yield_mm, 0),
deficit_35_var = ifelse(Yield_mm < thresh_35_var, thresh_35_var - Yield_mm, 0),
deficit_30_var = ifelse(Yield_mm < thresh_30_var, thresh_30_var - Yield_mm, 0),
deficit_25_var = ifelse(Yield_mm < thresh_25_var, thresh_25_var - Yield_mm, 0),
deficit_20_var = ifelse(Yield_mm < thresh_20_var, thresh_20_var - Yield_mm, 0),
deficit_15_var = ifelse(Yield_mm < thresh_15_var, thresh_15_var - Yield_mm, 0),
deficit_10_var = ifelse(Yield_mm < thresh_10_var, thresh_10_var - Yield_mm, 0),
deficit_05_var = ifelse(Yield_mm < thresh_05_var, thresh_05_var - Yield_mm, 0),
deficit_02_var = ifelse(Yield_mm < thresh_05_var, thresh_05_var - Yield_mm, 0))
# # fill missing dates
# dat_clean_join_deficit <- fill_missing_dates(dat_clean_join_deficit, dates = date, groups = site_name, pad_ends = TRUE)
#
# # fill ragged basin, subbasin, region, date variables
# dat_clean_join_deficit <- dat_clean_join_deficit %>%
# select(-c(basin, subbasin, region, CalendarYear, Month, MonthName, WaterYear, DayofYear, designation)) %>%
# left_join(siteinfo %>%
# mutate(designation = ifelse(site_name %in% unique(dat_clean_big$site_name), "big", "little")) %>%
# select(site_name, basin, subbasin, region, designation))
# dat_clean_join_deficit <- add_date_variables(dat_clean_join_deficit, dates = date, water_year_start = 10)
#
# # summarize by month
# defdur_month <- dat_clean_join_deficit %>%
# group_by(site_name, basin, subbasin, region, designation, CalendarYear, Month, MonthName, WaterYear) %>%
# summarize(ndays = n(),
# duration_50_fix = sum(deficit_50_fix > 0),
# duration_45_fix = sum(deficit_45_fix > 0),
# duration_40_fix = sum(deficit_40_fix > 0),
# duration_35_fix = sum(deficit_35_fix > 0),
# duration_30_fix = sum(deficit_30_fix > 0),
# duration_25_fix = sum(deficit_25_fix > 0),
# duration_20_fix = sum(deficit_20_fix > 0),
# duration_15_fix = sum(deficit_15_fix > 0),
# duration_10_fix = sum(deficit_10_fix > 0),
# duration_05_fix = sum(deficit_05_fix > 0),
# duration_02_fix = sum(deficit_02_fix > 0),
#
# duration_50_var = sum(deficit_50_var > 0),
# duration_45_var = sum(deficit_45_var > 0),
# duration_40_var = sum(deficit_40_var > 0),
# duration_35_var = sum(deficit_35_var > 0),
# duration_30_var = sum(deficit_30_var > 0),
# duration_25_var = sum(deficit_25_var > 0),
# duration_20_var = sum(deficit_20_var > 0),
# duration_15_var = sum(deficit_15_var > 0),
# duration_10_var = sum(deficit_10_var > 0),
# duration_05_var = sum(deficit_05_var > 0),
# duration_02_var = sum(deficit_02_var > 0),
#
# deficit_50_fix = sum(deficit_50_fix),
# deficit_45_fix = sum(deficit_45_fix),
# deficit_40_fix = sum(deficit_40_fix),
# deficit_35_fix = sum(deficit_35_fix),
# deficit_30_fix = sum(deficit_30_fix),
# deficit_25_fix = sum(deficit_25_fix),
# deficit_20_fix = sum(deficit_20_fix),
# deficit_15_fix = sum(deficit_15_fix),
# deficit_10_fix = sum(deficit_10_fix),
# deficit_05_fix = sum(deficit_05_fix),
# deficit_02_fix = sum(deficit_02_fix),
#
# deficit_50_var = sum(deficit_50_var),
# deficit_45_var = sum(deficit_45_var),
# deficit_40_var = sum(deficit_40_var),
# deficit_35_var = sum(deficit_35_var),
# deficit_30_var = sum(deficit_30_var),
# deficit_25_var = sum(deficit_25_var),
# deficit_20_var = sum(deficit_20_var),
# deficit_15_var = sum(deficit_15_var),
# deficit_10_var = sum(deficit_10_var),
# deficit_05_var = sum(deficit_05_var),
# deficit_02_var = sum(deficit_02_var)) %>%
# ungroup() %>%
# left_join(wateravail %>% select(basin, WaterYear, totalyield, totalyield_z))
# summarize by summer
defdur_ssn <- dat_clean_join_deficit %>%
filter(!is.na(Yield_mm), Month %in% c(7:9)) %>%
group_by(site_name, basin, subbasin, region, designation, CalendarYear, WaterYear) %>%
summarize(ndays = n(),
propdays = ndays/(31+31+30),
duration_50_fix_prop = sum(deficit_50_fix > 0) / ndays,
duration_45_fix_prop = sum(deficit_45_fix > 0) / ndays,
duration_40_fix_prop = sum(deficit_40_fix > 0) / ndays,
duration_35_fix_prop = sum(deficit_35_fix > 0) / ndays,
duration_30_fix_prop = sum(deficit_30_fix > 0) / ndays,
duration_25_fix_prop = sum(deficit_25_fix > 0) / ndays,
duration_20_fix_prop = sum(deficit_20_fix > 0) / ndays,
duration_15_fix_prop = sum(deficit_15_fix > 0) / ndays,
duration_10_fix_prop = sum(deficit_10_fix > 0) / ndays,
duration_05_fix_prop = sum(deficit_05_fix > 0) / ndays,
duration_02_fix_prop = sum(deficit_02_fix > 0) / ndays,
duration_50_var_prop = sum(deficit_50_var > 0) / ndays,
duration_45_var_prop = sum(deficit_45_var > 0) / ndays,
duration_40_var_prop = sum(deficit_40_var > 0) / ndays,
duration_35_var_prop = sum(deficit_35_var > 0) / ndays,
duration_30_var_prop = sum(deficit_30_var > 0) / ndays,
duration_25_var_prop = sum(deficit_25_var > 0) / ndays,
duration_20_var_prop = sum(deficit_20_var > 0) / ndays,
duration_15_var_prop = sum(deficit_15_var > 0) / ndays,
duration_10_var_prop = sum(deficit_10_var > 0) / ndays,
duration_05_var_prop = sum(deficit_05_var > 0) / ndays,
duration_02_var_prop = sum(deficit_02_var > 0) / ndays,
deficit_50_fix = sum(deficit_50_fix),
deficit_45_fix = sum(deficit_45_fix),
deficit_40_fix = sum(deficit_40_fix),
deficit_35_fix = sum(deficit_35_fix),
deficit_30_fix = sum(deficit_30_fix),
deficit_25_fix = sum(deficit_25_fix),
deficit_20_fix = sum(deficit_20_fix),
deficit_15_fix = sum(deficit_15_fix),
deficit_10_fix = sum(deficit_10_fix),
deficit_05_fix = sum(deficit_05_fix),
deficit_02_fix = sum(deficit_02_fix),
deficit_50_var = sum(deficit_50_var),
deficit_45_var = sum(deficit_45_var),
deficit_40_var = sum(deficit_40_var),
deficit_35_var = sum(deficit_35_var),
deficit_30_var = sum(deficit_30_var),
deficit_25_var = sum(deficit_25_var),
deficit_20_var = sum(deficit_20_var),
deficit_15_var = sum(deficit_15_var),
deficit_10_var = sum(deficit_10_var),
deficit_05_var = sum(deficit_05_var),
deficit_02_var = sum(deficit_02_var)) %>%
ungroup() %>%
filter(propdays >= 0.90) %>%
left_join(wateravail %>% select(basin, WaterYear, totalyield, totalyield_z))Create plotting functions
durationplotfun <- function(subbas, bigG, months, wateryears, dropsites = NA) {
# filter and summarize data
# dd <- defdur_month %>%
# filter(subbasin == subbas | site_name == bigG, Month %in% months, WaterYear %in% wateryears, !site_name %in% dropsites) %>%
# mutate(WaterYear = factor(WaterYear, levels = wateryears)) %>%
# group_by(site_name, designation, WaterYear, totalyield_z) %>%
# summarize(duration_50_fix_prop = sum(duration_50_fix) / sum(ndays),
# duration_45_fix_prop = sum(duration_45_fix) / sum(ndays),
# duration_40_fix_prop = sum(duration_40_fix) / sum(ndays),
# duration_35_fix_prop = sum(duration_35_fix) / sum(ndays),
# duration_30_fix_prop = sum(duration_30_fix) / sum(ndays),
# duration_25_fix_prop = sum(duration_25_fix) / sum(ndays),
# duration_20_fix_prop = sum(duration_20_fix) / sum(ndays),
# duration_15_fix_prop = sum(duration_15_fix) / sum(ndays),
# duration_10_fix_prop = sum(duration_10_fix) / sum(ndays),
# duration_05_fix_prop = sum(duration_05_fix) / sum(ndays),
# duration_02_fix_prop = sum(duration_02_fix) / sum(ndays)) %>%
# ungroup()
dd <- defdur_ssn %>%
filter(subbasin == subbas | site_name == bigG, WaterYear %in% wateryears, !site_name %in% dropsites) %>%
mutate(WaterYear = factor(WaterYear, levels = wateryears))
# order sites, Big G first
mysites <- c(unique(unlist(dd %>% filter(designation == "big") %>% select(site_name))),
unique(unlist(dd %>% filter(designation == "little") %>% select(site_name))))
# among site StDev ~ percentile
p_sds <- dd %>%
gather(duration_50_fix_prop:duration_02_fix_prop, key = "metric", value = "duration") %>%
mutate(quant = as.numeric(gsub(".*?([0-9]+).*", "\\1", metric)) ) %>%
filter(designation == "little") %>%
group_by(WaterYear, totalyield_z, metric, quant) %>%
summarize(sddur = sd(duration, na.rm = TRUE)) %>%
ungroup() %>%
left_join(dd %>% filter(site_name == bigG) %>% select(WaterYear, duration_50_fix_prop) %>% rename(dur50 = duration_50_fix_prop)) %>%
ggplot(aes(x = quant, y = sddur, color = dur50, group = WaterYear, shape = WaterYear)) +
geom_point() +
scale_color_gradient(low = "blue", high = "red", limits = c(0,1)) +
stat_smooth() +
xlab("Low flow threshold (percentile)") + ylab("Among-site SD(duration)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
legend.position = "bottom", legend.direction = "vertical", legend.key.height = unit(0.3, 'cm'))
# barplot 50th perc
p30 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = duration_50_fix_prop)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,1) +
ylab("Days below low flow threshold (%JAS)") + ggtitle("50th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 25th perc.
p20 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = duration_25_fix_prop)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,1) +
ylab("Days below low flow threshold (%JAS)") + ggtitle("25th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 10th perc.
p10 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = duration_10_fix_prop)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,1) +
ylab("Days below low flow threshold (%JAS)") + ggtitle("10th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 5th perc.
p05 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = duration_05_fix_prop)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,1) +
ylab("Days below low flow threshold (%JAS)") + ggtitle("5th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 2nd perc
p02 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = duration_02_fix_prop)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,1) +
ylab("Days below low flow threshold (%JAS)") + ggtitle("2nd perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# arrange plots
egg::ggarrange(p02 + theme(plot.margin = margin(r = 1, t = 5, b = 5)),
p05 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p10 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p20 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p30 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p_sds,
nrow = 1, widths = c(1,1,1,1,1,2.5))
}
deficitplotfun <- function(subbas, bigG, months, wateryears, dropsites = NA) {
# filter and summarize data
# dd <- defdur_month %>%
# filter(subbasin == subbas | site_name == bigG, Month %in% months, WaterYear %in% wateryears, !site_name %in% dropsites) %>%
# mutate(WaterYear = factor(WaterYear, levels = wateryears)) %>%
# group_by(site_name, designation, WaterYear, totalyield_z) %>%
# summarize(deficit_50_fix = sum(deficit_50_fix),
# deficit_45_fix = sum(deficit_45_fix),
# deficit_40_fix = sum(deficit_40_fix),
# deficit_35_fix = sum(deficit_35_fix),
# deficit_30_fix = sum(deficit_30_fix),
# deficit_25_fix = sum(deficit_25_fix),
# deficit_20_fix = sum(deficit_20_fix),
# deficit_15_fix = sum(deficit_15_fix),
# deficit_10_fix = sum(deficit_10_fix),
# deficit_05_fix = sum(deficit_05_fix),
# deficit_02_fix = sum(deficit_02_fix)) %>%
# ungroup()
dd_all <- defdur_ssn %>% filter(subbasin == subbas | site_name == bigG, WaterYear %in% wateryears, !site_name %in% dropsites)
dd <- defdur_ssn %>%
filter(subbasin == subbas | site_name == bigG, WaterYear %in% wateryears, !site_name %in% dropsites) %>%
mutate(WaterYear = factor(WaterYear, levels = wateryears))
# get y-axis limit
ymax <- max(dd %>% select(deficit_50_fix:deficit_02_fix))
# order sites, Big G first
mysites <- c(unique(unlist(dd %>% filter(designation == "big") %>% select(site_name))),
unique(unlist(dd %>% filter(designation == "little") %>% select(site_name))))
# among site StDev ~ percentile
p_sds <- dd %>%
gather(deficit_50_fix:deficit_02_fix, key = "metric", value = "deficit") %>%
mutate(quant = as.numeric(gsub(".*?([0-9]+).*", "\\1", metric)) ) %>%
filter(designation == "little") %>%
group_by(WaterYear, totalyield_z, metric, quant) %>%
summarize(sddur = sd(deficit, na.rm = TRUE)) %>%
ungroup() %>%
left_join(dd %>% filter(site_name == bigG) %>% select(WaterYear, deficit_50_fix) %>% rename(def50 = deficit_50_fix)) %>%
ggplot(aes(x = quant, y = sddur, color = def50, group = WaterYear, shape = WaterYear)) +
geom_point() +
scale_color_gradient(low = "blue", high = "red", limits = c(0,max(dd_all$deficit_50_fix))) +
stat_smooth() +
xlab("Low flow threshold (percentile)") + ylab("Among-site SD(deficit)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
legend.position = "bottom", legend.direction = "vertical", legend.key.height = unit(0.3, 'cm'))
# barplot 50th perc
p30 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = deficit_50_fix)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,ymax) +
ylab("Drought deficit (mm, JAS)") + ggtitle("50th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 25th perc.
p20 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = deficit_25_fix)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,ymax) +
ylab("Drought deficit (mm, JAS)") + ggtitle("25th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 10th perc.
p10 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = deficit_10_fix)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,ymax) +
ylab("Drought deficit (mm, JAS)") + ggtitle("10th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 5th perc.
p05 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = deficit_05_fix)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,ymax) +
ylab("Drought deficit (mm, JAS)") + ggtitle("5th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 2nd perc
p02 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = deficit_02_fix)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,ymax) +
ylab("Drought deficit (JAS)") + ggtitle("2nd perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# arrange plots
egg::ggarrange(p02 + theme(plot.margin = margin(r = 1, t = 5, b = 5)),
p05 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p10 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p20 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p30 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p_sds,
nrow = 1, widths = c(1,1,1,1,1,2.5))
}Show proportion of days (July - September) below different low flow thresholds (derived from long-term Big G data) for each site during a relatively wet year and a dry year. Then, for each year, plot the relationship between the among site (little g’s only) standard deviation of low flow duration and the low flow threshold used to calculate duration
durationplotfun(subbas = "West Brook", bigG = "South River Conway NWIS", months = c(7:9), wateryears = c(2020, 2021))
durationplotfun(subbas = "Staunton River", bigG = "Rapidan River NWIS", months = c(7:9), wateryears = c(2019, 2020))
durationplotfun(subbas = "Big Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9), wateryears = c(2020, 2021), dropsites = c("LangfordCreekUpper", "NicolaCreek", "WernerCreek", "LangfordCreekLower"))
durationplotfun(subbas = "Snake River", bigG = "Pacific Creek at Moran NWIS", months = c(7:9), wateryears = c(2021, 2022))
durationplotfun(subbas = "Shields River", bigG = "Yellowstone River Livingston NWIS", months = c(7:9), wateryears = c(2020, 2019), dropsites = c("Shields River Valley Ranch"))
durationplotfun(subbas = "Duck Creek", bigG = "Yellowstone River Livingston NWIS", months = c(7:9), wateryears = c(2021, 2019))
durationplotfun(subbas = "Donner Blitzen", bigG = "Donner Blitzen River nr Frenchglen NWIS", months = c(7:9), wateryears = c(2022, 2020))
Show total drought deficit (mm) relative to different low flow thresholds (derived from long-term Big G data) for each site during a relatively wet year and a dry year. Then, for each year, plot the relationship between the among site (little g’s only) standard deviation of deficit and the low flow threshold used to calculate deficit
deficitplotfun(subbas = "West Brook", bigG = "South River Conway NWIS", months = c(7:9), wateryears = c(2020, 2021))
deficitplotfun(subbas = "Staunton River", bigG = "Rapidan River NWIS", months = c(7:9), wateryears = c(2019, 2020))
deficitplotfun(subbas = "Big Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9), wateryears = c(2020, 2021), dropsites = c("LangfordCreekUpper", "NicolaCreek", "WernerCreek", "LangfordCreekLower"))
deficitplotfun(subbas = "Snake River", bigG = "Pacific Creek at Moran NWIS", months = c(7:9), wateryears = c(2021, 2022))
deficitplotfun(subbas = "Shields River", bigG = "Yellowstone River Livingston NWIS", months = c(7:9), wateryears = c(2020, 2019), dropsites = c("Shields River Valley Ranch"))
deficitplotfun(subbas = "Duck Creek", bigG = "Yellowstone River Livingston NWIS", months = c(7:9), wateryears = c(2021, 2019))
deficitplotfun(subbas = "Donner Blitzen", bigG = "Donner Blitzen River nr Frenchglen NWIS", months = c(7:9), wateryears = c(2022, 2020))
Organize data
dat_clean_sub_deficit <- dat_clean_sub %>%
mutate(deficit_50_fix = ifelse(z_Yield_mm < thresh_50_fix, abs(thresh_50_fix - z_Yield_mm), 0),
deficit_45_fix = ifelse(z_Yield_mm < thresh_45_fix, abs(thresh_45_fix - z_Yield_mm), 0),
deficit_40_fix = ifelse(z_Yield_mm < thresh_40_fix, abs(thresh_40_fix - z_Yield_mm), 0),
deficit_35_fix = ifelse(z_Yield_mm < thresh_35_fix, abs(thresh_35_fix - z_Yield_mm), 0),
deficit_30_fix = ifelse(z_Yield_mm < thresh_30_fix, abs(thresh_30_fix - z_Yield_mm), 0),
deficit_25_fix = ifelse(z_Yield_mm < thresh_25_fix, abs(thresh_25_fix - z_Yield_mm), 0),
deficit_20_fix = ifelse(z_Yield_mm < thresh_20_fix, abs(thresh_20_fix - z_Yield_mm), 0),
deficit_15_fix = ifelse(z_Yield_mm < thresh_15_fix, abs(thresh_15_fix - z_Yield_mm), 0),
deficit_10_fix = ifelse(z_Yield_mm < thresh_10_fix, abs(thresh_10_fix - z_Yield_mm), 0),
deficit_05_fix = ifelse(z_Yield_mm < thresh_05_fix, abs(thresh_05_fix - z_Yield_mm), 0),
deficit_02_fix = ifelse(z_Yield_mm < thresh_05_fix, abs(thresh_05_fix - z_Yield_mm), 0))
# summarize by summer
defdur_ssn_sub <- dat_clean_sub_deficit %>%
filter(!is.na(Yield_mm), Month %in% c(7:9)) %>%
group_by(site_name, basin, subbasin, region, designation, CalendarYear, WaterYear) %>%
summarize(ndays = n(),
propdays = ndays/(31+31+30),
duration_50_fix_prop = sum(deficit_50_fix > 0) / ndays,
duration_45_fix_prop = sum(deficit_45_fix > 0) / ndays,
duration_40_fix_prop = sum(deficit_40_fix > 0) / ndays,
duration_35_fix_prop = sum(deficit_35_fix > 0) / ndays,
duration_30_fix_prop = sum(deficit_30_fix > 0) / ndays,
duration_25_fix_prop = sum(deficit_25_fix > 0) / ndays,
duration_20_fix_prop = sum(deficit_20_fix > 0) / ndays,
duration_15_fix_prop = sum(deficit_15_fix > 0) / ndays,
duration_10_fix_prop = sum(deficit_10_fix > 0) / ndays,
duration_05_fix_prop = sum(deficit_05_fix > 0) / ndays,
duration_02_fix_prop = sum(deficit_02_fix > 0) / ndays,
deficit_50_fix = sum(deficit_50_fix),
deficit_45_fix = sum(deficit_45_fix),
deficit_40_fix = sum(deficit_40_fix),
deficit_35_fix = sum(deficit_35_fix),
deficit_30_fix = sum(deficit_30_fix),
deficit_25_fix = sum(deficit_25_fix),
deficit_20_fix = sum(deficit_20_fix),
deficit_15_fix = sum(deficit_15_fix),
deficit_10_fix = sum(deficit_10_fix),
deficit_05_fix = sum(deficit_05_fix),
deficit_02_fix = sum(deficit_02_fix)) %>%
ungroup() %>%
mutate(designation = ifelse(is.na(designation), "big", designation)) %>%
filter(propdays >= 0.90) %>%
left_join(wateravail %>% select(basin, WaterYear, totalyield, totalyield_z))Create plotting functions. These are the same as defined above, but instead grab the “defdur_ssn_sub” object for site-level low flow thresholds.
durationplotfun_sub <- function(subbas, bigG, months, wateryears, dropsites = NA) {
# filter and summarize data
dd <- defdur_ssn_sub %>%
filter(subbasin == subbas | site_name == bigG, WaterYear %in% wateryears, !site_name %in% dropsites) %>%
mutate(WaterYear = factor(WaterYear, levels = wateryears))
# order sites, Big G first
mysites <- c(unique(unlist(dd %>% filter(designation == "big") %>% select(site_name))),
unique(unlist(dd %>% filter(designation == "little") %>% select(site_name))))
# among site StDev ~ percentile
p_sds <- dd %>%
gather(duration_50_fix_prop:duration_02_fix_prop, key = "metric", value = "duration") %>%
mutate(quant = as.numeric(gsub(".*?([0-9]+).*", "\\1", metric)) ) %>%
filter(designation == "little") %>%
group_by(WaterYear, totalyield_z, metric, quant) %>%
summarize(sddur = sd(duration, na.rm = TRUE)) %>%
ungroup() %>%
left_join(dd %>% filter(site_name == bigG) %>% select(WaterYear, duration_25_fix_prop) %>% rename(dur25 = duration_25_fix_prop)) %>%
ggplot(aes(x = quant, y = sddur, color = dur25, group = WaterYear, shape = WaterYear)) +
stat_smooth() +
geom_point() +
scale_color_gradient(low = "blue", high = "red", limits = c(0,1)) +
xlab("Low flow threshold (percentile)") + ylab("Among-site SD(duration)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
legend.position = "bottom", legend.direction = "vertical", legend.key.height = unit(0.3, 'cm'))
# barplot 50th perc
p30 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = duration_50_fix_prop)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,1) +
ylab("Days below low flow threshold (%JAS)") + ggtitle("50th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 25th perc.
p20 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = duration_25_fix_prop)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,1) +
ylab("Days below low flow threshold (%JAS)") + ggtitle("25th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 10th perc.
p10 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = duration_10_fix_prop)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,1) +
ylab("Days below low flow threshold (%JAS)") + ggtitle("10th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 5th perc.
p05 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = duration_05_fix_prop)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,1) +
ylab("Days below low flow threshold (%JAS)") + ggtitle("5th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 2nd perc
p02 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = duration_02_fix_prop)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,1) +
ylab("Days below low flow threshold (%JAS)") + ggtitle("2nd perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# arrange plots
egg::ggarrange(p02 + theme(plot.margin = margin(r = 1, t = 5, b = 5)),
p05 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p10 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p20 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p30 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p_sds,
nrow = 1, widths = c(1,1,1,1,1,2.5))
}
deficitplotfun_sub <- function(subbas, bigG, months, wateryears, dropsites = NA) {
# filter and summarize data
dd_all <- defdur_ssn_sub %>% filter(subbasin == subbas | site_name == bigG, WaterYear %in% wateryears, !site_name %in% dropsites)
dd <- defdur_ssn_sub %>%
filter(subbasin == subbas | site_name == bigG, WaterYear %in% wateryears, !site_name %in% dropsites) %>%
mutate(WaterYear = factor(WaterYear, levels = wateryears))
# get y-axis limit
ymax <- max(dd %>% select(deficit_50_fix:deficit_02_fix))
# order sites, Big G first
mysites <- c(unique(unlist(dd %>% filter(designation == "big") %>% select(site_name))),
unique(unlist(dd %>% filter(designation == "little") %>% select(site_name))))
# among site StDev ~ percentile
p_sds <- dd %>%
gather(deficit_50_fix:deficit_02_fix, key = "metric", value = "deficit") %>%
mutate(quant = as.numeric(gsub(".*?([0-9]+).*", "\\1", metric)) ) %>%
filter(designation == "little") %>%
group_by(WaterYear, totalyield_z, metric, quant) %>%
summarize(sddur = sd(deficit, na.rm = TRUE)) %>%
ungroup() %>%
left_join(dd %>% filter(site_name == bigG) %>% select(WaterYear, deficit_25_fix) %>% rename(def25 = deficit_25_fix)) %>%
ggplot(aes(x = quant, y = sddur, color = def25, group = WaterYear, shape = WaterYear)) +
stat_smooth() +
geom_point() +
scale_color_gradient(low = "blue", high = "red", limits = c(0,max(dd_all$deficit_25_fix))) +
xlab("Low flow threshold (percentile)") + ylab("Among-site SD(deficit)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
legend.position = "bottom", legend.direction = "vertical", legend.key.height = unit(0.3, 'cm'))
# barplot 50th perc
p30 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = deficit_50_fix)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,ymax) +
ylab("Drought deficit (mm, JAS)") + ggtitle("50th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 25th perc.
p20 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = deficit_25_fix)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,ymax) +
ylab("Drought deficit (mm, JAS)") + ggtitle("25th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 10th perc.
p10 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = deficit_10_fix)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,ymax) +
ylab("Drought deficit (mm, JAS)") + ggtitle("10th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 5th perc.
p05 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = deficit_05_fix)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,ymax) +
ylab("Drought deficit (mm, JAS)") + ggtitle("5th perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# barplot 2nd perc
p02 <- dd %>%
ggplot(aes(x = factor(site_name, levels = (mysites)), y = deficit_02_fix)) +
geom_bar(aes(fill = designation), stat = "identity") +
scale_fill_manual(values = c("grey20", "grey55")) +
facet_wrap2(~WaterYear, ncol = 1, strip = strip_themed(background_x = elem_list_rect(fill = alpha(unique(layer_data(p_sds)[,1]), 0.5)))) +
ylim(0,ymax) +
ylab("Drought deficit (JAS)") + ggtitle("2nd perc.") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
legend.position = "none", axis.title.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# arrange plots
egg::ggarrange(p02 + theme(plot.margin = margin(r = 1, t = 5, b = 5)),
p05 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p10 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p20 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p30 + theme(axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = margin(r = 1, l = 1)),
p_sds,
nrow = 1, widths = c(1,1,1,1,1,2.5))
}Show proportion of days (July - September) below different low flow thresholds (derived from temporally restricted, site-specific data) for each site during a relatively wet year and a dry year. Then, for each year, plot the relationship between the among site (little g’s only) standard deviation of low flow duration and the low flow threshold used to calculate duration
durationplotfun_sub(subbas = "West Brook", bigG = "South River Conway NWIS", months = c(7:9), wateryears = c(2022, 2020, 2021))
durationplotfun_sub(subbas = "Staunton River", bigG = "Rapidan River NWIS", months = c(7:9), wateryears = c(2019, 2021, 2020))
durationplotfun_sub(subbas = "Big Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9), wateryears = c(2020, 2019))
durationplotfun_sub(subbas = "Coal Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9), wateryears = c(2019, 2021))
durationplotfun_sub(subbas = "Snake River", bigG = "Pacific Creek at Moran NWIS", months = c(7:9), wateryears = c(2021, 2020, 2022))
durationplotfun_sub(subbas = "Shields River", bigG = "Yellowstone River Livingston NWIS", months = c(7:9), wateryears = c(2020, 2019, 2023), dropsites = c("Shields River Valley Ranch"))
durationplotfun_sub(subbas = "Duck Creek", bigG = "Yellowstone River Livingston NWIS", months = c(7:9), wateryears = c(2021, 2022, 2017))
durationplotfun_sub(subbas = "Donner Blitzen", bigG = "Donner Blitzen River nr Frenchglen NWIS", months = c(7:9), wateryears = c(2022, 2020))
Show total drought deficit (mm) relative to different low flow thresholds (derived from long-term Big G data) for each site during a relatively wet year and a dry year. Then, for each year, plot the relationship between the among site (little g’s only) standard deviation of deficit and the low flow threshold used to calculate deficit
deficitplotfun_sub(subbas = "West Brook", bigG = "South River Conway NWIS", months = c(7:9), wateryears = c(2022, 2020, 2021))
deficitplotfun_sub(subbas = "Staunton River", bigG = "Rapidan River NWIS", months = c(7:9), wateryears = c(2019, 2021, 2020))
deficitplotfun_sub(subbas = "Big Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9), wateryears = c(2020, 2019))
deficitplotfun_sub(subbas = "Coal Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9), wateryears = c(2019, 2021))
deficitplotfun_sub(subbas = "Snake River", bigG = "Pacific Creek at Moran NWIS", months = c(7:9), wateryears = c(2021, 2020, 2022))
deficitplotfun_sub(subbas = "Shields River", bigG = "Yellowstone River Livingston NWIS", months = c(7:9), wateryears = c(2020, 2019, 2023), dropsites = c("Shields River Valley Ranch"))
deficitplotfun_sub(subbas = "Duck Creek", bigG = "Yellowstone River Livingston NWIS", months = c(7:9), wateryears = c(2021, 2022, 2017))
deficitplotfun_sub(subbas = "Donner Blitzen", bigG = "Donner Blitzen River nr Frenchglen NWIS", months = c(7:9), wateryears = c(2022, 2020))
Purpose: Show what high-resolution temporal data (hourly) reveals about network diversity in streamflow response to individual storms (peak flow magnitude and timing and recession rates) and during low flow (diel fluctuations)
The Wedge Model for the West Brook
Purpose: At coarser time scales (summarized by event/baseflow periods), show how streamflow heterogeneity expands and contracts during wet to dry periods.
Purpose: Show what high-resolution spatial data reveals about network diversity in streamflow at a single point in time.
Purpose: Explore the effect of groundwater on relative summer (July-September) water availability.
Load PASTA daily derived parameters: summarize as July-September site-specific means (across all years)
pasta <- read_csv("C:/Users/jbaldock/OneDrive - DOI/Documents/USGS/EcoDrought/EcoDrought Working/EcoDrought-Analysis/Covariates/pasta_derived_parameters_daily.csv") %>%
mutate(Month = month(date)) %>%
rename(CalendarYear = year) %>%
filter(Month %in% c(7:9)) %>%
group_by(site_name) %>%
summarize(meanRatio = mean(meanRatio, na.rm = TRUE),
phaseLag = mean(phaseLag, na.rm = TRUE),
amplitudeRatio = mean(amplitudeRatio, na.rm = TRUE))
pasta# A tibble: 73 × 4
site_name meanRatio phaseLag amplitudeRatio
<chr> <dbl> <dbl> <dbl>
1 Avery Brook 0.891 2.31 0.264
2 BigCreekLower 0.936 2.22 0.376
3 BigCreekMiddle 0.840 1.29 0.411
4 BigCreekUpper 0.664 2.31 0.179
5 Buck Creek 0.699 1.19 0.353
6 CoalCreekHeadwaters 0.670 2.11 0.209
7 CoalCreekLower 0.843 2.91 0.503
8 CoalCreekMiddle 0.692 1.27 0.223
9 CoalCreekNorth 0.733 2.49 0.246
10 Crandall Creek 0.758 1.76 0.612
# ℹ 63 more rowsCreate flow by groundwater plotting function.
mdaystib <- tibble(Month = c(1:12), mdays = c(31,28,31,30,31,30,31,31,30,31,30,31))
gwflowfun <- function (subbas, years, dropsites, months = c(1:12)) {
dat_clean %>%
filter(subbasin == subbas, CalendarYear %in% years, Month %in% months) %>%
group_by(site_name, subbasin, designation, CalendarYear) %>% #, Month, MonthName) %>%
summarise(ss = n(),
logYield = mean(logYield, na.rm = TRUE)) %>%
ungroup() %>%
#left_join(mdaystib) %>%
mutate(pdays = ss/92#,
#YearMonth = paste(CalendarYear, "_", Month, sep = "")
) %>%
filter(pdays > 0.9,
!site_name %in% dropsites) %>%
group_by(CalendarYear) %>%
mutate(z_logYield = scale(logYield, center = TRUE, scale = TRUE)[,1]) %>%
ungroup() %>%
left_join(pasta) %>%
ggplot(aes(x = amplitudeRatio, y = z_logYield)) +
geom_abline(intercept = 0, slope = 0, linetype = 2) +
geom_smooth(method = "lm", color = "black") +
geom_point(aes(color = site_name)) +
facet_wrap(~CalendarYear, nrow = 1) +
#facet_wrap2(~CalendarYear, nrow = 1, ncol = 5, trim_blank = FALSE) +
#facet_grid(cols = vars(Month), rows = vars(CalendarYear)) +
theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
}Plot the relationship between standardized annual summer mean discharge and amplitude ratio from PASTA, where lower amplitude ratio values are indicative of greater groundwater availability. Mean flow for each site is standardized by year to remove interannual variation in climate/regional water availability
gwflowfun(subbas = "West Brook", dropsites = c("West Brook Reservoir", "Mitchell Brook"), years = c(2020:2024))
gwflowfun(subbas = "Staunton River", dropsites = NA, years = c(2019:2022))
gwflowfun(subbas = "Snake River", dropsites = NA, years = c(2018, 2020:2022), months = c(7:9))
gwflowfun(subbas = "Shields River", dropsites = NA, years = c(2017,2019,2020,2022,2023), months = c(7:9))
Create violin plot function. Don’t love this plotting style, but shows that data reasonably well
defdurplotfun <- function(subbas, bigG, months) {
dd_little <- defdur_month %>% filter(subbasin == subbas, Month %in% months, designation == "little")
myyrs <- unique(dd_little$WaterYear)
dd_little <- dd_little %>% mutate(WaterYear = factor(WaterYear, levels = (sort(myyrs))))
dd_big <- defdur_month %>%
filter(site_name == bigG, Month %in% months) %>%
mutate(WaterYear = factor(WaterYear, levels = (sort(myyrs)))) %>%
filter(WaterYear %in% unique(myyrs))
m <- 7
p_def_7 <- ggplot() +
geom_violin(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, fill = totalyield_z), trim = TRUE, scale = "width") +
scale_fill_gradient2(low = "red", mid = "white", high = "blue", midpoint = 0) +
geom_line(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, group = site_name, color = site_name), orientation = "y") +
geom_point(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, group = site_name, color = site_name)) +
geom_line(data = dd_big %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, group = site_name), orientation = "y") +
geom_point(data = dd_big %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, group = site_name), size = 2) +
coord_flip() +
theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p_dur_7 <- ggplot() +
geom_violin(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, fill = totalyield_z), trim = TRUE, scale = "width") +
scale_fill_gradient2(low = "red", mid = "white", high = "blue", midpoint = 0) +
geom_line(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, group = site_name, color = site_name), orientation = "y") +
geom_point(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, group = site_name, color = site_name)) +
geom_line(data = dd_big %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, group = site_name), orientation = "y") +
geom_point(data = dd_big %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, group = site_name), size = 2) +
coord_flip() +
theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
m <- 8
p_def_8 <- ggplot() +
geom_violin(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, fill = totalyield_z), trim = TRUE, scale = "width") +
scale_fill_gradient2(low = "red", mid = "white", high = "blue", midpoint = 0) +
geom_line(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, group = site_name, color = site_name), orientation = "y") +
geom_point(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, group = site_name, color = site_name)) +
geom_line(data = dd_big %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, group = site_name), orientation = "y") +
geom_point(data = dd_big %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, group = site_name), size = 2) +
coord_flip() +
theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p_dur_8 <- ggplot() +
geom_violin(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, fill = totalyield_z), trim = TRUE, scale = "width") +
scale_fill_gradient2(low = "red", mid = "white", high = "blue", midpoint = 0) +
geom_line(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, group = site_name, color = site_name), orientation = "y") +
geom_point(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, group = site_name, color = site_name)) +
geom_line(data = dd_big %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, group = site_name), orientation = "y") +
geom_point(data = dd_big %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, group = site_name), size = 2) +
coord_flip() +
theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
m <- 9
p_def_9 <- ggplot() +
geom_violin(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, fill = totalyield_z), trim = TRUE, scale = "width") +
scale_fill_gradient2(low = "red", mid = "white", high = "blue", midpoint = 0) +
geom_line(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, group = site_name, color = site_name), orientation = "y") +
geom_point(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, group = site_name, color = site_name)) +
geom_line(data = dd_big %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, group = site_name), orientation = "y") +
geom_point(data = dd_big %>% filter(Month == m), aes(y = WaterYear, x = deficit_10_fix, group = site_name), size = 2) +
coord_flip() +
theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p_dur_9 <- ggplot() +
geom_violin(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, fill = totalyield_z), trim = TRUE, scale = "width") +
scale_fill_gradient2(low = "red", mid = "white", high = "blue", midpoint = 0) +
geom_line(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, group = site_name, color = site_name), orientation = "y") +
geom_point(data = dd_little %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, group = site_name, color = site_name)) +
geom_line(data = dd_big %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, group = site_name), orientation = "y") +
geom_point(data = dd_big %>% filter(Month == m), aes(y = WaterYear, x = duration_10_fix, group = site_name), size = 2) +
coord_flip() +
theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
return(ggarrange(p_def_7, p_def_8, p_def_9,
p_dur_7, p_dur_8, p_dur_9,
nrow = 2, ncol = 3, common.legend = TRUE, legend = "right") %>%
annotate_figure(top = "July August September"))
}defdurplotfun(subbas = "West Brook", bigG = "South River Conway NWIS", months = c(7:9))defdurplotfun(subbas = "Paine Run", bigG = "South River Harriston NWIS", months = c(7:9))defdurplotfun(subbas = "Staunton River", bigG = "Rapidan River NWIS", months = c(7:9))defdurplotfun(subbas = "Big Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9))defdurplotfun(subbas = "Coal Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9))defdurplotfun(subbas = "McGee Creek", bigG = "North Fork Flathead River NWIS", months = c(7:9))defdurplotfun(subbas = "Snake River", bigG = "Pacific Creek at Moran NWIS", months = c(7:9))defdurplotfun(subbas = "Shields River", bigG = "Yellowstone River Livingston NWIS", months = c(7:9))defdurplotfun(subbas = "Duck Creek", bigG = "Yellowstone River Livingston NWIS", months = c(7:9))defdurplotfun(subbas = "Donner Blitzen", bigG = "Donner Blitzen River nr Frenchglen NWIS", months = c(7:9))How does the spatial variation in stream flow compare to temporal variation in streamflow at a single reference location?
Create plotting function
myplotfun_spacetime <- function(subbas) {
# format data
data_sub <- dat_clean %>% filter(subbasin == subbas)
sitesdesig <- data_sub %>% group_by(site_name) %>% summarize(designation = unique(designation), ss = n()) %>% ungroup()
data_sub_com <- data_sub %>%
select(date, site_name, logYield) %>%
spread(key = site_name, value = logYield) %>%
drop_na() %>%
gather(key = site_name, value = logYield, 2:ncol(.))
data_sub_com_little <- data_sub_com %>% filter(site_name %in% unlist(sitesdesig %>% filter(designation == "little") %>% select(site_name)))
data_sub_com_big <- data_sub_com %>% filter(site_name %in% unlist(sitesdesig %>% filter(designation == "big") %>% select(site_name)))
# calculate spatial and temporal variation
temporal_sd <- sd(data_sub_com_big$logYield)
relsd_ts <- data_sub_com_little %>%
group_by(date) %>%
summarize(little_sdlf = sd(logYield)) %>%
ungroup() %>%
mutate(big_temporal_sdlf = temporal_sd,
relsdlf = little_sdlf - big_temporal_sdlf) %>%
left_join(data_sub_com_big)
# plot time series
par(mfrow = c(1,2), mar = c(3,3,0.5,0.5), mgp = c(2,0.6,0))
plot(relsdlf ~ date, relsd_ts, xlab = "Date", ylab = "Spatial SD - temporal SD")
abline(h = 0, col = "red", lty = 2)
#par(mar = c(3,3,0.1,0.1), mgp = c(2.5,0.6,0))
legend("topright", legend = subbas, bty = "n")
# plot sd by big G flow
plot(relsdlf ~ logYield, relsd_ts, axes = FALSE, xlab = "log(Yield at big G)", ylab = "")
axis(1)
axis(2, labels = NA)
box(bty = "O")
abline(h = 0, col = "red", lty = 2)
}myplotfun_spacetime(subbas = "West Brook")myplotfun_spacetime(subbas = "Paine Run")myplotfun_spacetime(subbas = "Staunton River")myplotfun_spacetime(subbas = "Big Creek")myplotfun_spacetime(subbas = "Coal Creek")myplotfun_spacetime(subbas = "McGee Creek")myplotfun_spacetime(subbas = "Snake River")myplotfun_spacetime(subbas = "Shields River")myplotfun_spacetime(subbas = "Duck Creek")myplotfun_spacetime(subbas = "Donner Blitzen")Load drought/climate metrics data and view temporal range of 1-month SPI and DSCI.
drought <- read_csv("C:/Users/jbaldock/OneDrive - DOI/Documents/USGS/EcoDrought/EcoDrought Working/EcoDrought-Analysis/Explore Data/EcoD_Basin_MonthlyDroughtMetrics.csv") %>% rename(Month = month, CalendarYear = year) #%>% select(-date)
datatable(drought %>% select(basin, date, spi_1, dsci_monthly_1) %>% gather(key = "metric", value = "value", spi_1:dsci_monthly_1) %>% drop_na() %>% group_by(basin, metric) %>% summarize(mindate = min(date), maxdate = max(date)))Show distribtuions of SPI and DSCI during recent years, by basin
# SPI
drought %>%
filter(CalendarYear >= 2016) %>%
select(basin, CalendarYear, Month, date, spi_1, spi_6, spi_12, spi_24) %>%
gather(key = "metric", value = "value", spi_1:spi_24) %>%
mutate(metric = factor(metric, levels = c("spi_1", "spi_6", "spi_12", "spi_24"))) %>%
ggplot(aes(x = basin, y = value, fill = metric)) +
geom_boxplot() +
xlab("") + ylab("Standardized Precipitation Index (2016-2023)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
# DSCI
drought %>%
filter(CalendarYear >= 2016) %>%
select(basin, CalendarYear, Month, date, dsci_monthly_1, dsci_monthly_6, dsci_monthly_12, dsci_monthly_24) %>%
gather(key = "metric", value = "value", dsci_monthly_1:dsci_monthly_24) %>%
mutate(metric = factor(metric, levels = c("dsci_monthly_1", "dsci_monthly_6", "dsci_monthly_12", "dsci_monthly_24"))) %>%
ggplot(aes(x = basin, y = value, fill = metric)) +
geom_boxplot() +
xlab("") + ylab("Drought Severity and Coverage Index (2016-2024)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
Create plotting function
myplotfun_mmq <- function(subbas) {
# summarize monthly mean yield by site, year
tt <- dat_clean %>%
mutate(basin = ifelse(subbasin == "Duck Creek", "Shields River", basin)) %>%
filter(!is.na(logYield)) %>%
group_by(region, basin, subbasin, site_name, designation, CalendarYear, Month) %>%
summarize(meanQ = mean(logYield, na.rm = TRUE), nobs = n()) %>%
ungroup() %>%
filter(nobs >= 27) %>%
left_join(drought)
# filter to focal subbasin, keep only sites/months with at least 3 years of data
tt2 <- tt %>% filter(subbasin == subbas) %>% mutate(keep = paste(site_name, Month, sep = "_"))
keeps <- tt2 %>%
group_by(site_name, Month) %>%
summarize(nobs = n()) %>%
filter(nobs >= 3) %>%
mutate(keep = paste(site_name, Month, sep = "_"))
tt2_little <- tt2 %>% filter(keep %in% keeps$keep, designation == "little", Month %in% c(6:10))
tt2_big <- tt2 %>% filter(keep %in% keeps$keep, designation == "big", Month %in% c(6:10))
# plot SPI-9
p1 <- ggplot() +
geom_point(data = tt2_little, aes(x = spi_9, y = meanQ, color = site_name, group = site_name)) +
geom_smooth(data = tt2_little, aes(x = spi_9, y = meanQ, color = site_name, group = site_name), method = "lm", se = FALSE) +
geom_point(data = tt2_big, aes(x = spi_9, y = meanQ), color = "black") +
geom_smooth(data = tt2_big, aes(x = spi_9, y = meanQ), color = "black", method = "lm", se = FALSE) +
facet_wrap(~Month, nrow = 1) +
theme_bw() + theme(legend.position = "none", panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
xlab("Standardized precipitation index, 9 month") + ylab("Mean monthly log(yield)")
# plot DSCI-1
p2 <- ggplot() +
geom_point(data = tt2_little, aes(x = dsci_monthly_1, y = meanQ, color = site_name, group = site_name)) +
geom_smooth(data = tt2_little, aes(x = dsci_monthly_1, y = meanQ, color = site_name, group = site_name), method = "lm", se = FALSE) +
geom_point(data = tt2_big, aes(x = dsci_monthly_1, y = meanQ), color = "black") +
geom_smooth(data = tt2_big, aes(x = dsci_monthly_1, y = meanQ), color = "black", method = "lm", se = FALSE) +
facet_wrap(~Month, nrow = 1) +
theme_bw() + theme(legend.position = "none", panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
xlab("Drought severity and coverage index, 1 month") + ylab("Mean monthly log(yield)")
# combine plots
ggarrange(p1, p2, nrow = 2)
}myplotfun_mmq(subbas = "West Brook")
myplotfun_mmq(subbas = "Paine Run")
myplotfun_mmq(subbas = "Staunton River")
myplotfun_mmq(subbas = "Big Creek")
myplotfun_mmq(subbas = "Coal Creek")
myplotfun_mmq(subbas = "McGee Creek")
myplotfun_mmq(subbas = "Snake River")
myplotfun_mmq(subbas = "Shields River")
myplotfun_mmq(subbas = "Duck Creek")
myplotfun_mmq(subbas = "Donner Blitzen")
myplotfun_dist <- function(subbas, distance = c("euclidean", "manhattan")) {
# filter by subbasin
data_sub <- dat_clean %>%
mutate(basin = ifelse(subbasin == "Duck Creek", "Shields River", basin)) %>%
filter(subbasin == subbas)
# little g's
data_sub_little <- data_sub %>%
filter(designation == "little") %>%
select(site_name, basin, subbasin, date, CalendarYear, Month, logYield)
# big G
data_sub_big <- data_sub %>%
filter(designation == "big") %>%
select(basin, subbasin, date, CalendarYear, Month, logYield) %>%
rename(logYield_big = logYield)
# join little and big
data_sub_join <- data_sub_little %>% full_join(data_sub_big)
# calculate time series similarity metrics
dist_sum <- data_sub_join %>%
group_by(basin, site_name, CalendarYear, Month) %>%
summarize(numdays = n(),
dist_ts = ifelse(distance == "euclidean",
EuclideanDistance(logYield, logYield_big),
ManhattanDistance(logYield, logYield_big))) %>%
ungroup() %>%
left_join(drought) %>%
filter(numdays >= 27, Month %in% c(6:10)) %>%
mutate(keep = paste(site_name, Month, sep = "_"))
# filter to focal subbasin, keep only sites/months with at least 3 years of data
keeps <- dist_sum %>%
group_by(site_name, Month) %>%
summarize(nobs = n()) %>%
filter(nobs >= 3) %>%
mutate(keep = paste(site_name, Month, sep = "_"))
dist_sum <- dist_sum %>% filter(keep %in% keeps$keep)
# Plots
p1 <- dist_sum %>%
ggplot(aes(x = spi_9, y = dist_ts, color = site_name, group = site_name)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
facet_wrap(~Month, nrow = 1) +
theme_bw() + theme(legend.position = "none", panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
xlab("Standardized precipitation index, 9 month") + ylab("Distance")
p2 <- dist_sum %>%
ggplot(aes(x = dsci_monthly_1, y = dist_ts, color = site_name, group = site_name)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
facet_wrap(~Month, nrow = 1) +
theme_bw() + theme(legend.position = "none", panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
xlab("Drought severity and coverage index, 1 month") + ylab("Distance")
# combine plots
ggarrange(p1, p2, nrow = 2)
}myplotfun_dist(subbas = "West Brook", distance = "euclidean")Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
myplotfun_dist(subbas = "Paine Run", distance = "euclidean")Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
myplotfun_dist(subbas = "Staunton River", distance = "euclidean")Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
myplotfun_dist(subbas = "Big Creek", distance = "euclidean")Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
myplotfun_dist(subbas = "Coal Creek", distance = "euclidean")Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
myplotfun_dist(subbas = "McGee Creek", distance = "euclidean")Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
myplotfun_dist(subbas = "Snake River", distance = "euclidean")Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
myplotfun_dist(subbas = "Shields River", distance = "euclidean")Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
myplotfun_dist(subbas = "Duck Creek", distance = "euclidean")Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
myplotfun_dist(subbas = "Donner Blitzen", distance = "euclidean")Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
Error : There are missing values in the series.
# summarize monthly mean yield by site, year
tt <- dat_clean %>%
mutate(basin = ifelse(subbasin == "Duck Creek", "Shields River", basin)) %>%
filter(!is.na(logYield), designation == "little") %>%
group_by(region, basin, subbasin, site_name, designation, CalendarYear, Month) %>%
summarize(meanQ = mean(logYield, na.rm = TRUE), nobs = n()) %>%
ungroup() %>%
filter(nobs >= 27) %>%
left_join(drought)
# calculate correlations
mylist_dsci_slope <- list()
mylist_dsci_pval <- list()
mylist_dsci_rho <- list()
mylist_spi_slope <- list()
mylist_spi_pval <- list()
mylist_spi_rho <- list()
# list of sites to iterate over
sites <- unique(tt$site_name)
for(i in 1:length(sites)) {
d <- tt %>% filter(site_name == sites[i])
mymat_dsci_slope <- matrix(data = NA, nrow = 12, ncol = 24)
mymat_dsci_pval <- matrix(data = NA, nrow = 12, ncol = 24)
mymat_dsci_rho <- matrix(data = NA, nrow = 12, ncol = 24)
mymat_spi_slope <- matrix(data = NA, nrow = 12, ncol = 24)
mymat_spi_pval <- matrix(data = NA, nrow = 12, ncol = 24)
mymat_spi_rho <- matrix(data = NA, nrow = 12, ncol = 24)
for(j in 1:12) {
for(k in 1:24) {
if(nrow(d %>% filter(Month == j)) < 3) next
# DSCI
d2 <- d %>% filter(Month == j) %>% select(meanQ, paste("dsci_monthly_", k, sep = ""))
## slope and p-value
try(mylm <- lm(unlist(d2[,1]) ~ unlist(d2[,2])), silent = TRUE)
try(mymat_dsci_slope[j,k] <- summary(mylm)$coefficients[2,1], silent = TRUE)
try(mymat_dsci_pval[j,k] <- summary(mylm)$coefficients[2,4], silent = TRUE)
## pearons r correlation
try(mymat_dsci_rho[j,k] <- cor(d2[,1], d2[,2], method = "pearson", use = "complete.obs"), silent = TRUE)
# SPI
d2 <- d %>% filter(Month == j) %>% select(meanQ, paste("spi_", k, sep = ""))
## slope and p-value
try(mylm <- lm(unlist(d2[,1]) ~ unlist(d2[,2])), silent = TRUE)
try(mymat_spi_slope[j,k] <- summary(mylm)$coefficients[2,1], silent = TRUE)
try(mymat_spi_pval[j,k] <- summary(mylm)$coefficients[2,4], silent = TRUE)
## pearons r correlation
try(mymat_spi_rho[j,k] <- cor(d2[,1], d2[,2], method = "pearson", use = "complete.obs"), silent = TRUE)
}
}
mylist_dsci_slope[[i]] <- mymat_dsci_slope
mylist_dsci_pval[[i]] <- mymat_dsci_pval
mylist_dsci_rho[[i]] <- mymat_dsci_rho
mylist_spi_slope[[i]] <- mymat_spi_slope
mylist_spi_pval[[i]] <- mymat_spi_pval
mylist_spi_rho[[i]] <- mymat_spi_rho
print(i)
}
# generate plots/heatmaps of correlation
n <- 100
for(i in 1:length(sites)) {
# streamflow mean ~ DSCI
jpeg(paste("./Climate Sensitivity/Climate sensitivity plots/Streamflow_DSCI_correlation_TimeVarying_", "site", i, "_mean.jpg", sep = ""), height = 6, width = 7, units = "in", res = 500)
filled.contour(x = 1:24, y = 1:12, z = t(corlist_dsci_mean[[i]]),
levels = seq(from = -1, to = 1, length = n+1), col = hcl.colors(n, "PRGn"),
plot.title = title(main = sites[i],
xlab = "DSCI - accumulation time scale (months)",
ylab = "Month"),
plot.axes = {
axis(1, at = seq(from = 1, to = 23, by = 2))
axis(2, at = 1:12, labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"))
contour(x = 1:24, y = 1:12, z = t(corlist_dsci_mean[[i]]), add = TRUE, labcex = 1, levels = seq(from = -1, to = 1, length = 11))
})
dev.off()
# streamflow mean ~ SPI
jpeg(paste("./Climate Sensitivity/Climate sensitivity plots/Streamflow_SPI_correlation_TimeVarying_", "site", i, "_mean.jpg", sep = ""), height = 6, width = 7, units = "in", res = 500)
filled.contour(x = 1:24, y = 1:12, z = t(corlist_spi_mean[[i]]),
levels = seq(from = -1, to = 1, length = n+1), col = hcl.colors(n, "PRGn"),
plot.title = title(main = sites[i],
xlab = "SPI - accumulation time scale (months)",
ylab = "Month"),
plot.axes = {
axis(1, at = seq(from = 1, to = 23, by = 2))
axis(2, at = 1:12, labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"))
contour(x = 1:24, y = 1:12, z = t(corlist_spi_mean[[i]]), add = TRUE, labcex = 1, levels = seq(from = -1, to = 1, length = 11))
})
dev.off()
# streamflow minimum ~ DSCI
jpeg(paste("./Climate Sensitivity/Climate sensitivity plots/Streamflow_DSCI_correlation_TimeVarying_", "site", i, "_min.jpg", sep = ""), height = 6, width = 7, units = "in", res = 500)
filled.contour(x = 1:24, y = 1:12, z = t(corlist_dsci_min[[i]]),
levels = seq(from = -1, to = 1, length = n+1), col = hcl.colors(n, "PRGn"),
plot.title = title(main = sites[i],
xlab = "DSCI - accumulation time scale (months)",
ylab = "Month"),
plot.axes = {
axis(1, at = seq(from = 1, to = 23, by = 2))
axis(2, at = 1:12, labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"))
contour(x = 1:24, y = 1:12, z = t(corlist_dsci_min[[i]]), add = TRUE, labcex = 1, levels = seq(from = -1, to = 1, length = 11))
})
dev.off()
# streamflow minimum ~ SPI
jpeg(paste("./Climate Sensitivity/Climate sensitivity plots/Streamflow_SPI_correlation_TimeVarying_", "site", i, "_min.jpg", sep = ""), height = 6, width = 7, units = "in", res = 500)
filled.contour(x = 1:24, y = 1:12, z = t(corlist_spi_min[[i]]),
levels = seq(from = -1, to = 1, length = n+1), col = hcl.colors(n, "PRGn"),
plot.title = title(main = sites[i],
xlab = "SPI - accumulation time scale (months)",
ylab = "Month"),
plot.axes = {
axis(1, at = seq(from = 1, to = 23, by = 2))
axis(2, at = 1:12, labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"))
contour(x = 1:24, y = 1:12, z = t(corlist_spi_min[[i]]), add = TRUE, labcex = 1, levels = seq(from = -1, to = 1, length = 11))
})
dev.off()
}
# plot raw data
plot(meanQ ~ spi_9, tt %>% filter(site_name == "Big Creek NWIS", Month == 10))
plot(z_flow_mean_monthly ~ spi_7, dat3 %>% filter(site_name == "Donner Blitzen River nr Frenchglen NWIS", month == 6))
filled.contour(x = 1:24, y = 1:12, z = t(mylist_dsci_slope[[1]]),
levels = seq(from = -1, to = 1, length = 50),
col = hcl.colors(50, "PRGn"),
plot.title = title(main = sites[i],
xlab = "DSCI - accumulation time scale (months)",
ylab = "Month"),
plot.axes = {
axis(1, at = seq(from = 1, to = 23, by = 2))
axis(2, at = 1:12, labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"))
contour(x = 1:24, y = 1:12, z = t(mylist_dsci_slope[[i]]), add = TRUE, labcex = 1, levels = seq(from = -1, to = 1, length = 11))
})
filled.contour(x = 1:24, y = 1:12, z = t(mylist_spi_pval[[1]]),
levels = seq(from = 0, to = 1, length = 100),
col = hcl.colors(100, "Blues"))
filled.contour(x = 1:24, y = 1:12, z = t(mylist_spi_slope[[1]]),
levels = seq(from = -1, to = 1, length = 100),
col = hcl.colors(100, "PRGn"))